{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting the necessary data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You just need to download this ~28 MB file only once"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-08-26 16:50:33--  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265.filt.fastq.gz\n",
      "           => ‘SRR003265.filt.fastq.gz’\n",
      "Resolving ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)... 193.62.192.8\n",
      "Connecting to ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)|193.62.192.8|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /vol1/ftp/phase3/data/NA18489/sequence_read ... done.\n",
      "==> SIZE SRR003265.filt.fastq.gz ... 28919712\n",
      "==> PASV ... done.    ==> RETR SRR003265.filt.fastq.gz ... done.\n",
      "Length: 28919712 (28M) (unauthoritative)\n",
      "\n",
      "SRR003265.filt.fast 100%[===================>]  27.58M  6.22MB/s    in 4.4s    \n",
      "\n",
      "2018-08-26 16:50:39 (6.22 MB/s) - ‘SRR003265.filt.fastq.gz’ saved [28919712]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!rm -f SRR003265.filt.fastq.gz 2>/dev/null\n",
    "!wget -nd ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265.filt.fastq.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The recipe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tiago/anaconda3/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n",
      "/home/tiago/anaconda3/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n"
     ]
    }
   ],
   "source": [
    "from collections import defaultdict\n",
    "import gzip\n",
    "\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SRR003265.31 SRR003265.31 3042NAAXX:3:1:1252:1819 length=51 GGGAAAAGAAAAACAAACAAACAAAAACAAAACACAGAAACAAAAAAACCA\n",
      "{'phred_quality': [40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 30, 23, 40, 32, 35, 29, 40, 16, 40, 40, 32, 35, 31, 40, 40, 39, 22, 40, 24, 20, 28, 31, 12, 31, 10, 22, 28, 13, 26, 20, 23, 23]}\n"
     ]
    }
   ],
   "source": [
    "recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')\n",
    "rec = next(recs)\n",
    "print(rec.id, rec.description, rec.seq)\n",
    "print(rec.letter_annotations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A: 28.60 7411965\n",
      "T: 29.58 7666885\n",
      "G: 20.68 5359334\n",
      "N: 0.14 37289\n",
      "C: 21.00 5444053\n"
     ]
    }
   ],
   "source": [
    "recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')\n",
    "cnt = defaultdict(int)\n",
    "for rec in recs:\n",
    "    for letter in rec.seq:\n",
    "        cnt[letter] += 1\n",
    "tot = sum(cnt.values())\n",
    "for letter, cnt in cnt.items():\n",
    "    print('%s: %.2f %d' % (letter, 100. * cnt / tot, cnt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 51)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAIMCAYAAAD1pfEjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3WtwW+l95/nfA4DgBQApAgRJtdTduhFU2912bLf7Lqovie3Ozoyz2WTLu5PE5bjKVXF2N5mdrd1k3qQ2M9lkXux4NjWVTGU2TjxJZhNXJjPJZm0nnpZaVLfV3e5uX1rdIkFJLbUuBAnwCt5AXJ59AUCmJEriBcDBOfh+qlQiDw8O/rS7JP70PM//b6y1AgAAAACgWfmcLgAAAAAAgLshuAIAAAAAmhrBFQAAAADQ1AiuAAAAAICmRnAFAAAAADQ1gisAAAAAoKkRXAEAAAAATY3gCgAAAABoagRXAAAAAEBTI7gCAAAAAJpawOkC7qavr88eOHDA6TIAAAAAAHXw1ltvZay18Xvd19TB9cCBA3rzzTedLgMAAAAAUAfGmMtbuY+twgAAAACApkZwBQAAAAA0NYIrAAAAAKCpEVwBAAAAAE2N4AoAAAAAaGoEVwAAAABAUyO4AgAAAACaGsEVAAAAANDUCK4AAAAAgKZGcAUAAAAANDWCKwAAAACgqRFcAQAAAABNjeAKAAAAAGhqBFcAAAAAQFMjuAIAAAAAmhrBFQAAAADQ1AiuAAAAAICmRnAFAAAAADQ1gisAAAAAoKkRXAEAAByQKxT12X/zir51NuV0KQDQ9AiuAAAADnjr8px+cHVBf/nWVadLAYCmR3AFAABwwOmJjCTpzIWM1gslh6sBgOZGcAUAAHDAaDKtjjaflteLevuDOafLAYCmRnAFAABosMxSTu9eX9TnnzyggM9oNJl2uiQAaGoEVwAAgAZ79Xx5m/BPPrJXH3+gV6MTBFcAuBuCKwAAQIONJjPa09Wmh/f1aCTRp7PXFpXO5pwuCwCaFsEVAACggay1Oj2R1tNH+uT3GY0k4pKkV86z6goAd0JwBQAAaKDk1JKmszmNDPVJkh6+r0fRUFCjyYzDlQFA8yK4AgAANNDpynnWZ4bKK60+n9EzR/p0eiKtUsk6WRoANC2CKwAAQAONTmR0OB7Svj2dN66NJOLKLK3rvclFBysDgOZFcAUAAGiQtXxRr1+c0bHKamtVddsw3YUBYHMEVwAAgAZ589KccoWSRhJ9N13v7+7Q0cEI81wB4A4IrgAAAA1yeiKtNr/R4wdjt33teCKuty7PaTlXcKAyAGhuBFcAAIAGGZ3I6NEHowq1B2772kgirnzR6syFGQcqA4DmRnAFAABogOnsms5NLurYLduEqx490KvONj/nXAFgEwRXAACABnj1fHlO68gtjZmq2gN+PXk4xjlXANgEwRUAAKABTiczioaC+tDe7jveMzLUp0szK7o8s9zAygCg+RFcAQAA6sxaq9GJjJ450iefz9zxvpFEeTWWVVcAuBnBFQAAoM7GUllllnI6NrT5+daqg30h7e/t1KlkpkGVAYA7EFwBAADq7HSl4dKxO5xvrTLGaCQR15kLGa0XSo0oDQBcgeAKAABQZ6cnMkoMhDXY03HPe0eG4lpeL+rtD+YaUBkAuAPBFQAAoI7W8kW9/v7sPVdbq546EpPfZzjnCgAbEFwBAADq6I33Z7VeKN3zfGtVd0ebPv7AHua5AsAGBFcAAIA6Oj2RVtDv0+MHY1t+zchQXGevLSqzlKtjZQDgHgRXAACAOjo9kdEnD/aqM+jf8muqY3FemaC7MABIBFcAAIC6mV5c01gqu+XzrVWP7OtRNBTknCsAVBBcAQAA6uR0ZcV0q+dbq3w+o2eO9Gl0IqNSydajNABwFYIrAABAnYxOpNUXDuqhwe5tv3YkEVdmKadzqcU6VAYA7kJwBQAAqINSyeqViYyeOdInn89s+/UjlVXaU2wXBgCCKwAAQD28N7momeX1bZ9vrerv7tDRwQjnXAFABFcAAIC62On51o2OJ+J66/KclnOFWpUFAK5EcAUAAKiD0xNpHR2MqL+7Y8fPGEnElS9anbkwU8PKAMB9CK4AAAA1trJe0JuX5na12ipJjx7oVWebX6MTbBcG0NoIrgAAADX2+vuzWi+Wdny+tao94NcTh6KccwXQ8giuAAAANXY6mVF7wKfHDkZ3/ayRRFyXZlb0wcxKDSoDAHciuAIAANTY6Ym0HjsYVUebf9fPOp4or9qeYrswgBZGcAUAAKihyYVVTUwvaWSX24SrDvaFtL+3k+3CAFoawRUAAKCGbozBSeyuMVOVMUYjibjOXJhRvliqyTMBwG0IrgAAADV0eiKjeKRdwwORmj1zZCiupVxBb1+eq9kzAcBNCK4AAAA1UipZvTKR1rGhPhljavbcp47E5PcZnWK7MIAWRXAFAACokXevL2puJV+z861V3R1t+vgDe5jnCqBlEVwBAABqpBosnz5Sm/OtG40MxXX22qIyS7maPxsAmh3BFQAAoEZOT6T1ob3dikfaa/7skcpYnFcqzZ8AoJUQXAEAQMsoFEt168y7nCvorctzNesmfKuH9/Wot6uNsTgAWhLBFQAAtIz/4T98T//o37yqxbV8zZ/9+vszyhdtzc+3Vvl9Rs8MxTU6kVGpZOvyHgDQrAiuAACgZbx5eU7nJhf1S3/6ltYLtV15HU1m1NHm0yce7K3pczcaGepTZimnc6nFur0HADQjgisAAGgJM0s5ZZZyeuxAVK+en9Gv/dUPZW3tVi5PT6T1+MGYOtr8NXvmrY5XzrmOJjnnCqC1EFwBAEBLGJ/KSpL+xxeO6J/8eEJ/9fY1feW/TNTk2dfmV3UhvaxjQ/U531rV392ho4MRzrkCaDkBpwsAAABohGSqHFyHByJ65kifrs6t6HdfmtD+PZ36bz95/66e/UplDE618289HU/E9dVX39dyrqBQOz/KAWgNrLgCAICWMD61pD1dbYpH2mWM0f/x04/o2FCffv0/vaNTu1zBHJ3IaKC7XUP94RpVe2cjibjyRavXLs7U/b0AoFkQXAEAQEtITmWVGIjIGCNJavP79Hv/+OMa6g/ry3/6lt69vrCj5xZLVq+ez+jYUPzGs+vp0QO96mzzs10YQEshuAIAAM+z1iqZyuroYOSm65GONv3xFx5Td2ebfvGPv6vr86vbfvbZawuaX8nX/XxrVXvArycORXe9SgwAbkJwBQAAnnd9YU3ZXEGJgchtXxvs6dAffeGTWskV9YU/+u62Z7xWVz6fOdKY4CqVtwtfmlnRBzMrDXtPAHASwRUAAHjejcZMg7cHV0k6Otit3/+5T+hCeklf/tO3tzXj9fRERg/v61Ys3F6TWrei2gTq1ASrrgBaA8EVAAB4XnUUTqJ/8+AqSc8M9el3/puP6JXzGf36X72zpRmv2bW83v5gTseG6t9NeKNDfSHt29PJOVcALYMe6gAAwPOSqawGuzvU09V21/t+5hP7dW1uVV/5L0nt7+3UP/mJxF3vf+3irAolq5EGB1djjEYScf2/P7iufLGkNj9rEQC8jT/lAACA542lsnfcJnyr/+mFI/rZT+zX//XShL7+5pW73nt6Iq2uoF8ff3BPLcrcluOJPi3lCnr78lzD3xsAGo3gCgAAPK1QLOl8emnLwXXjjNd/9lfv3HU77umJjJ44FFN7wF+rcrfsqSN98vuMRjnnCqAFEFwBAICnXZ5d0XqhtGlH4Tupzng90h/Wl//sbb13ffG2e67Mruj9zHLDxuDcqrujTR9/YI9GkxlH3h8AGmnLwdUY4zfGfM8Y87eVzw8aY143xkwYY/7CGBOsXG+vfH6+8vUDG57x65Xr48aYT9f6mwEAALjVjY7C2wiuUnnG6x994ZMKtwf0i3/8XU0u3Dzj9fREOTA2ujHTRiNDcZ29vqCZpZxjNQBAI2xnxfVXJJ3b8Pm/lPQVa+2QpDlJX6xc/6KkOWvtEUlfqdwnY8yHJH1O0oclfUbS7xljGr+vBgAAtJTxqayMkY70h7f92r09nfqjL3xSS7nCbTNeT0+kdV9Phw7HQ7Usd1tGEnFZK71ynlVXAN62peBqjNkv6b+S9H9XPjeSnpf0l5Vbvibppyoff7byuSpff6Fy/2cl/bm1NmetfV/SeUmP1eKbAAAAuJPxVFYHYiF1Bnf27+UP7e3W7//cx3V++kczXgvFkl49n9GxobjKP+Y44+F9PertatMpxuIA8Litrrj+a0n/q6TqNO6YpHlrbaHy+VVJ+yof75N0RZIqX1+o3H/j+iavAQAAqIvxqawSA9tfbd3o2FBcv/3Tj9yY8fqDqwtaXCvoWMKZ861Vfp/RM0NxjSYzKpXuPXcWANzqnsHVGPMPJE1ba9/aeHmTW+09vna312x8vy8ZY940xryZTvOvhwAAYOfW8kVdyixv+3zrZn720fv1qz8+pP/49lX9069/X8ZITx92NrhK0shQnzJLOZ1L3d5ACgC8Yisrrk9L+kfGmEuS/lzlLcL/WtIeY0ygcs9+SdcrH1+VdL8kVb7eI2l24/VNXnODtfYPrLWPWmsfjceda3YAAADc70J6SSUrJbY4CudefuWFIf3MJ/br0syKPrKvR72hYE2euxsfe6BXkpScyjpcCQDUzz2Dq7X21621+621B1RurnTCWvuPJZ2U9DOV2z4v6a8rH/9N5XNVvn7CWmsr1z9X6Tp8UNKQpDdq9p0AAADcohrmarHiKpVnvP72Tz+iX3jyQf3Ss4dr8szd6guXw/Pscv4edwKAewXufcsd/W+S/twY8y8kfU/SH1au/6GkPzHGnFd5pfVzkmStfdcY83VJ70kqSPpla21xF+8PAABwV2OprIJ+nw701a7zb5vfp9/87MM1e95udXe0ye8zml1mJA4A79pWcLXWvizp5crHF7VJV2Br7Zqkn73D639L0m9tt0gAAICdSKayOhQPqc2/nQmA7uLzGfV2BVlxBeBp3v1THAAAtLzk1JKGa3S+tZlFQ22suALwNIIrAADwpOxaXtfmV5Wo0fnWZhYNBTW7vO50GQBQNwRXAADgScmpJUm1a8zUzGKhdoIrAE8juAIAAE8aT1U6CrfAVuHeUBvBFYCnEVwBAIAnJaeyCgX92ren0+lS6i4aatf8al7FknW6FACoC4IrAADwpPFUVkMDEfl8xulS6i4WCspaaX6FVVcA3kRwBQAAnpScyrbE+VZJ6g0FJYntwgA8i+AKAAA8J7OU08zyuhItcL5VKq+4StIMwRWARxFcAQCA51QbMx1tkeAarQTXOYIrAI8iuAIAAM+pBtdWmOEq/Si4suIKwKsIrgAAwHOSU1lFQ0H1hYNOl9IQvV2ccQXgbQRXAADgOeNTWSUGwjLG+x2FJSkY8CnSHiC4AvAsgisAAPAUa62SqayODnY7XUpDRcNBgisAzyK4AgAAT7k6t6rl9WLLnG+tioaCmmOOKwCPIrgCAABPSU6VGzMND4YdrqSxol1BzSwRXAF4E8EVAAB4yngluA614IorW4UBeBXBFQAAeEoyldV9PR3q7mhzupSGioaDml1Zl7XW6VIAoOYIrgAAwFPGp5Y0PNhaq61SeavweqGk5fWi06UAQM0RXAEAgGfkiyVdmF5SohWDa6gyy5VzrgA8iOAKAAA84/LMstaLJQ232PlWSYqFK8GVzsIAPIjgCgAAPGM8tSRJLTcKR5J6uyrBdTnncCUAUHsEVwAA4BnjU1n5jHSkv7VG4UhSLNQuSYzEAeBJBFcAAOAZyVRWB/pC6mjzO11Kw0UrW4Xn2CoMwIMIrgAAwDPGp7Iteb5VkkJBv4J+n2aY5QrAgwiuAADAE9byRV2aWW7J862SZIxRNBSkqzAATyK4AgAATzg/vSRr1ZIzXKuioSBbhQF4EsEVAAB4wngqK6k1OwpXRUNBtgoD8CSCKwAA8ITkVFbBgE8HYl1Ol+KYaCioWYIrAA8iuAIAAE8YS2V1JB5WwN+6P94QXAF4Vev+yQ4AADwlOZVt6fOtUjm4ZtcKWi+UnC4FAGqK4AoAAFxvYTWvyYW1lj7fKpWDq8QsVwDeQ3AFAACuNzFVbsw0PBh2uBJnVYMr24UBeA3BFQAAuN74jeDa7XAlziK4AvAqgisAAHC98VRW4faA7uvpcLoUR8UqwZWROAC8huAKAABcbzyVVWIgLGOM06U4qrd6xpXgCsBjCK4AAMDVrLV0FK7o7QrKGFZcAXgPwRUAALhaeimnuZV8y3cUliS/z2hPZ5tml3NOlwIANUVwBQAArpZMLUkSK64VvaGg5pbzTpcBADVFcAUAAK42llqUJA2z4iqp3KBphhVXAB5DcAUAAK6WnMqqLxxULNzudClNIcqKKwAPIrgCAABXG59a4nzrBtFQkOZMADyH4AoAAFyrVLKaoKPwTaKhoOZW1lUqWadLAYCaIbgCAADXuja/qpX1IudbN4iG2lUsWWXXCk6XAgA1Q3AFAACuNZbKSpISrLjeEA21SRINmgB4CsEVAAC4VnKqHFyH+sMOV9I8oqFyk6pZzrkC8BCCKwAAcK3xVFb79nQq0tHmdClNIxYKSiK4AvAWgisAAHCt5FRWR9kmfJNegisADyK4AgAAV8oXS7qQXuJ86y2qK66MxAHgJQRXAADgSu9nlpUvWjoK36Kjza+uoF9zBFcAHkJwBQAArjRe7ShMcL1Nb1eQrcIAPIXgCgAAXCk5lZXfZ3S4P+R0KU0nFg6yVRiApxBcAQCAK42nsjrYF1J7wO90KU2ntyuouRWCKwDvILgCAABXSk5lOd96B7FQUDNLBFcA3kFwBQAArrOyXtDl2RXOt95BNMQZVwDeQnAFAACuc356SdZKw4Nhp0tpSr2hoFbzRa2uF50uBQBqguAKAABch47Cd1ed5TrLOVcAHkFwBQAArpOcyqo94NODMToKbyZaDa6ccwXgEQRXAADgOuNTSxoaCMvvM06X0pSirLgC8BiCKwAAcJ3x1CLbhO/iRnBdzjlcCQDUBsEVAAC4yvzKuqYWc4zCuYtYqF2SGIkDwDMIrgAAwFWSU0uSpMQgwfVOIh0B+X1Gc2wVBuARBFcAAOAq41PljsJHCa535PMZ9XYxyxWAdxBcAQCAqyRTWUU6Ahrs7nC6lKYWCwXZKgzAMwiuAADAVcZTWQ0PRGQMHYXvpjfUxlZhAJ5BcAUAAK5hrdX4VJbzrVsQC7Vrhq3CADyC4AoAAFxjOpvTwmqejsJbEA1xxhWAdxBcAQCAa4ynyo2ZhllxvafeUFALq3kViiWnSwGAXSO4AgAA10hWOgonWHG9p1goKGul+dW806UAwK4RXAEAgGuMp7KKR9oVDQWdLqXpVf83mmO7MAAPILgCAADXGJ/Kcr51i6rBlQZNALyA4AoAAFyhVLJKTmXZJrxF1eBKgyYAXkBwBQAArnBlbkVr+ZKO0phpSwiuALyE4AoAAFyh2lGYGa5b09tFcAXgHQRXAADgCtWOwkP9YYcrcYdgwKdIR4DgCsATCK4AAMAVxlJZ3R/tVKg94HQprhENBQmuADyB4AoAAFwhOZXV8EC302W4CsEVgFcQXAEAQNNbL5R0Mb2s4UG2CW9HLBRkHA4ATyC4AgCApvd+ZlmFkmUUzjb1dgU1R3AF4AEEVwAA0PTGK42ZhukovC3RcHmrsLXW6VIAYFcIrgAAoOmNpxYV8Bkd6mOr8HbEQkGtF0tayhWcLgUAdoXgCgAAmt7rF2f10N5uBQP86LId1Vmuc8t5hysBgN3hT38AANDU5pbX9fYHc3ruaL/TpbhOLFwOrjPLOYcrAYDdIbgCAICmdiqZVslKzxNcty0aapckRuIAcD2CKwAAaGonxqbVFw7qI/t6nC7FdaKVrcIEVwBuR3AFAABNq1As6eXxaT073C+fzzhdjutEwwRXAN5AcAUAAE3r7Q/mtbhWYJvwDoWCfgUDPoIrANe7Z3A1xnQYY94wxvzAGPOuMeZ/r1w/aIx53RgzYYz5C2NMsHK9vfL5+crXD2x41q9Xro8bYz5dr28KAAB4w4mxaQV8RseG+pwuxZWMMYp2BQmuAFxvKyuuOUnPW2s/KunHJH3GGPOEpH8p6SvW2iFJc5K+WLn/i5LmrLVHJH2lcp+MMR+S9DlJH5b0GUm/Z4zx1/KbAQAA3nJibEqPHYwq0tHmdCmuFQ0RXAG43z2Dqy1bqnzaVvllJT0v6S8r178m6acqH3+28rkqX3/BGGMq1//cWpuz1r4v6bykx2ryXQAAAM+5Orei5NQS24R3KRYOaobgCsDltnTG1RjjN8Z8X9K0pG9LuiBp3lpbqNxyVdK+ysf7JF2RpMrXFyTFNl7f5DUAAAA3OTk2LYkxOLvV2xXU3ArBFYC7bSm4WmuL1tofk7Rf5VXShza7rfL7Zi3/7F2u38QY8yVjzJvGmDfT6fRWygMAAB700ti0DsS6dCgedroUV4uGgppdIrgCcLdtdRW21s5LelnSE5L2GGMClS/tl3S98vFVSfdLUuXrPZJmN17f5DUb3+MPrLWPWmsfjcfj2ykPAAB4xOp6UWcuzOg5Vlt3LRoKKpsrKFcoOl0KAOzYVroKx40xeyofd0r6cUnnJJ2U9DOV2z4v6a8rH/9N5XNVvn7CWmsr1z9X6Tp8UNKQpDdq9Y0AAADv+M6FjHKFkl44OuB0Ka4XDZVnuc6v5B2uBAB2LnDvW7RX0tcqHYB9kr5urf1bY8x7kv7cGPMvJH1P0h9W7v9DSX9ijDmv8krr5yTJWvuuMebrkt6TVJD0y9Za/ukPAADc5qWxaYWCfj12MOp0Ka4XqwTXmaV1DXR3OFwNAOzMPYOrtfaHkj62yfWL2qQrsLV2TdLP3uFZvyXpt7ZfJgAAaBXWWp0cm9YzQ30KBrZ1qgmb6K0EV0biAHAz/jYAAABN5dxkVpMLa2wTrpHqiussnYUBuBjBFQAANJWT4+UxOM8epUljLVTPuM4u5RyuBAB2juAKAACayomxaT2yr0f9Ec5j1sKerqCMkWZpzgTAxQiuAACgacwur+vtD+b0PGNwasbvM9rT2abZZVZcAbgXwRUAADSNU8lpWSuCa41FQ0GaMwFwNYIrAABoGifG0uoLt+uRfT1Ol+IpBFcAbkdwBQAATaFQLOnU+LSeG47L5zNOl+MpBFcAbkdwBQAATeGty3NaXCuwTbgOoqF2gisAVyO4AgCApnBifFptfqNnhvqcLsVzoqE2za3kVSpZp0sBgB0huAIAgKZw4ty0HjsYVaSjzelSPCcaalexZLW4xkgcAO5EcAUAAI67MruiieklPTfMNuF6iIWCkqQZtgsDcCmCKwAAcNzJ8WlJjMGpl95KcJ0juAJwKYIrAABw3EvnpnWwL6RD8bDTpXgSK64A3I7gCgAAHLWyXtCZizNsE66jaCW40lkYgFsRXAEAgKO+c35G64US24TriOAKwO0IrgAAwFEvjU0rFPTrsYNRp0vxrI42v7qCfoIrANciuAIAAMdYa/Xy+LSODcUVDPBjST31dgUJrgBci78hAACAY85NZjW5sMY24QaIhQmuANyL4AoAABxzYmxKkvTs0bjDlXhfNERwBeBeBFcAAOCYE2PT+sj+HvVHOpwuxfOibBUG4GIEVwAA4IjZ5XV978o8Y3AahBVXAG5GcAUAAI54eXxa1kovPERwbYRoOKjVfFGr60WnSwGAbSO4AgAAR5wYm1ZfuF0P39fjdCktIdpVnuU6s5xzuBIA2D6CKwAAaLh8saTRZFrPDcfl8xmny2kJ0VA5uM4t5x2uBAC2j+AKAAAa7q3Lc1pcK7BNuIFiYVZcAbgXwRUAADTcybFptfmNnhliDE6j9Fa2CtOgCYAbEVwBAEDDnRib1mMHowq3B5wupWXEQu2SCK4A3IngCgAAGurK7Iomppf0/NEBp0tpKd2dAfl9huAKwJUIrgAAoKFOjE1Lkp4/yvnWRjLGqLcrqLkVgisA9yG4AgCAhnppbFoH+0I62BdyupSWEwsFNbNEcAXgPgRXAADQMCvrBb12cYbVVodEQ0G2CgNwJYIrAABomFfPz2i9UCK4OiQaCmqWrcIAXIjgCgAAGubE2JTC7QF98kDU6VJaEiuuANyK4AoAABrCWquTY2kdG+pTMMCPIE6IhoKaX8mrUCw5XQoAbAt/awAAgIZ4b3JRqcU1Pcc2YcdEQ0FJ0vxq3uFKAGB7CK4AAKAhTpwrj8F5djjucCWtqxpc2S4MwG0IrgAAoCFOjE/ro/t71B/pcLqUllUNrozEAeA2BFcAAFB3M0s5ff/KPNuEHVYNrnN0FgbgMgRXAABQdy+Pp2WtGIPjsFh1xZWtwgBchuAKAADq7sT4tOKRdj18X4/TpbS0PV2VM65sFQbgMgRXAABQV/liSaPJtJ4bjsvnM06X09KCAZ8iHQG2CgNwHYIrAACoqzcvzSm7VmCbcJOIhYJsFQbgOgRXAABQVyfHp9XmN3pmiDE4zaA3FNTscs7pMgBgWwiuAACgrk6MTevxgzGF2wNOlwKVV1xnl/NOlwEA20JwBQAAdfPBzIrOTy8xBqeJRFlxBeBCBFcAAFA3J8enJTEGp5mUtwqvy1rrdCkAsGUEVwAAUDenkmk9GOvSwb6Q06WgIhYKKl+0WsoVnC4FALaM4AoAAOoiVyjqzIUZjdCUqalEQ+2SpFk6CwNwEYIrAACoi7cvz2s1X9SxoT6nS8EG0VCbJDESB4CrEFwBAEBdjE6kFfAZPXk45nQp2KC64jpHcAXgIgRXAABQF6cn0vr4A72KdLQ5XQo2iIWCklhxBeAuBFcAAFBzmaWczl5b1EiCbcLNprcSXDnjCsBNCK4AAKDmXj2fkSQdozFT0wkF/QoGfGwVBuAqBFcAAFBzp5Jp7elq08P7epwuBbcwxigWCrJVGICrEFwBAEBNWWt1eiKjZ470ye8zTpeDTfR2BdkqDMBVCK4AAKCmxlJZpbM55rc2sViY4ArAXQiuAACgpk5PpCVJx2jM1LRYcQXgNgRXAABQU6PJjIb6w9rb0+l0KbiDaChIcyYArkJwBQAANbO6XtQbl2Y1kmCbcDOLhYLK5grKFYpOlwIAW0JwBQAANfPGpVlm34fgAAAgAElEQVStF0o6NsQ24WZWneU6t5x3uBIA2BqCKwAAqJnRZFrBgE+PH4w5XQruIlYJrpxzBeAWBFcAAFAzpyfSeuxAVJ1Bv9Ol4C6iBFcALkNwBQAANTG5sKrk1JJG6Cbc9KrBdWY553AlALA1BFcAAFATpycykqRjzG9tetEbZ1xZcQXgDgRXAABQE6PJtOKRdh0djDhdCu5hT1dQxrBVGIB7EFwBAMCuFUtWr5zP6NhQn4wxTpeDe/D7jPZ0tmmG4ArAJQiuAABg185eW9D8Sl7Hmd/qGtFQUHMrBFcA7kBwBQAAu3Z6Ii1JevoIjZncIhZq18wSwRWAOxBcAQDAro0mM/rwfd3qC7c7XQq2qDfUxhlXAK5BcAUAALuSXcvr7Q/mNMI2YVeJhtrZKgzANQiuAABgV167OKtCyerYENuE3SQWCmpuJa9SyTpdCgDcE8EVAADsymgyra6gX594sNfpUrANvaGgiiWrhdW806UAwD0RXAEAwK6cnkjriUMxtQf8TpeCbYiFgpKkWbYLA3ABgisAANixD2ZWdGlmRSNsE3ad3mpwpUETABcguAIAgB0brYzBOUZjJteprrgyEgeAGxBcAQDAjo0m09q3p1OH+kJOl4JtilaCK52FAbgBwRUAAOxIvljSmQszGkn0yRjjdDnYpihbhQG4CMEVAADsyPevzCubK2hkiG3CbtTR5ldX0M9WYQCuQHAFAAA7cjqZls9ITx2mMZNbRUNBtgoDcAWCKwAA2JFTExn92P171NPV5nQp2KFoKKgZtgoDcAGCKwAA2Lb5lXX98Oq8jrFN2NWioaBml3NOlwEA90RwBQAA2/bq+RlZK40wBsfVoqGg5pbzTpcBAPdEcAUAANs2mkwr0hHQR/f3OF0KdiHaFdQMK64AXIDgCgAAtsVaq9MTaT1zpE8BPz9KuFk0HNRavqTV9aLTpQDAXfG3DQAA2JYL6SVdX1jjfKsHxCqzXFl1BdDsCK4AAGBbRpMZSdKxIcbguF1vVzm4ztJZGECTI7gCAIBtGZ1I61BfSPdHu5wuBbsUCxNcAbjDPYOrMeZ+Y8xJY8w5Y8y7xphfqVyPGmO+bYyZqPzeW7lujDG/a4w5b4z5oTHm4xue9fnK/RPGmM/X79sCAAD1kCsU9drFGVZbPSIaapdEcAXQ/Lay4lqQ9E+ttQ9JekLSLxtjPiTp1yS9ZK0dkvRS5XNJelHSUOXXlyT9vlQOupJ+Q9Ljkh6T9BvVsAsAANzhzUtzWsuXGIPjEVG2CgNwiXsGV2vtpLX27crHWUnnJO2T9FlJX6vc9jVJP1X5+LOS/r0te03SHmPMXkmflvRta+2stXZO0rclfaam3w0AAKir0Ym02vxGTxyKOV0KaqC7M6CAzxBcATS9bZ1xNcYckPQxSa9LGrDWTkrlcCupv3LbPklXNrzsauXana4DAACXGE1m9IkHexVqDzhdCmrAGKPeUJDgCqDpbTm4GmPCkv6jpF+11i7e7dZNrtm7XL/1fb5kjHnTGPNmOp3eankAAKDOprNrOje5yBgcj4l2BTVDcAXQ5LYUXI0xbSqH1j+z1v5V5fJUZQuwKr9PV65flXT/hpfvl3T9LtdvYq39A2vto9baR+Nx/mIEAKBZvHq+PAbnOOdbPSUaCmqO4AqgyW2lq7CR9IeSzllr/9WGL/2NpGpn4M9L+usN13+h0l34CUkLla3EfyfpU8aY3kpTpk9VrgEAABcYTWYUCwX1ob3dTpeCGoqyVRiAC2zlgMrTkn5e0jvGmO9Xrv0zSb8j6evGmC9K+kDSz1a+9g1JPynpvKQVSV+QJGvtrDHmn0v6buW+37TWztbkuwAAAHVVKlmdnsjomaE++Xybnf6BW0VDbBUG0PzuGVytta9o8/OpkvTCJvdbSb98h2d9VdJXt1MgAABw3rnUojJLOc63elA0FNTCal6FYkkB/7b6dgJAw/CnEwAAuKfTE+XzrSNDfQ5XglqLhsqzXOdW8g5XAgB3RnAFAAD3NJpM6+hgRP3dHU6XghqrBlfOuQJoZgRXAABwVyvrBb15aU4jdBP2pBjBFYALEFwBAMBdvX5xVuvFko6xTdiTegmuAFyA4AoAAO5qdCKt9oBPnzwQdboU1MGPVlxzDlcCAHdGcAUAAHc1mkzr8UMxdbT5nS4FdfCjFVeaMwFoXgRXAABwR9fmV3UhvUw3YQ9r8/sU6Qiw4gqgqRFcAQDAHb0ykZYkGjN5XCwU1AxnXAE0MYIrAAC4o9FkRoPdHRrqDztdCuooGgpqboXgCqB5EVwBAMCmiiWrV85ndGyoT8YYp8tBHUVDQc0sEVwBNC+CKwAA2NQPr85rYTWvY2wT9rxoKMg4HABNjeAKAAA2dXoiI2OkZ47QmMnroqF2za2sy1rrdCkAsCmCKwAA2NRoMq1H9vUoWhmXAu+KhtqUL1plcwWnSwGATRFcAQDAbRbX8vrelXmNDLFNuBVEQ+2SpFnOuQJoUgRXAABwm7cuzalYsnqabcItIVZZVZ+lszCAJkVwBQAAt7kytyJJOtwfcrgSNEJvNbiy4gqgSRFcAQDAbVILawr4jPoqW0jhbay4Amh2BFcAAHCb1MKaBro75PMxv7UV3FhxZSQOgCZFcAUAALdJLa5poJvV1lYRCvoVDPgIrgCaFsEVAADcJrW4pr09nU6XgQYxxigWChJcATQtgisAALiJtfbGVmG0jt4ugiuA5kVwBQAAN8nmClpZL2pvD8G1lcTCQc0QXAE0KYIrAAC4SWphTZI0QHBtKdFQUHMEVwBNiuAKAABuUg2urLi2FrYKA2hmBFcAAHCTanAd5IxrS4mFglrKFZQrFJ0uBQBuQ3AFAAA3SS2Wg2s/43BaSjRcnuU6t5x3uBIAuB3BFQAA3GRyYU2xUFDtAb/TpaCBBiLlFfYPZlccrqTs9Ysz+g+vf+B0GQCaBMEVAADcZGpxTYOcb205nzwQlc9IpyfSTpciSfo//z6p3/zbd1UqWadLAdAECK4AAOAmkwtrnG9tQT1dbfrYA706lXQ+uC6s5vXWB3Nay5d0Za45VoABOIvgCgAAbsKKa+t6NhHXD68uKLOUc7SOVyYyKlZWWsdTWUdrAdAcCK4AAOCGtXxRs8vrrLi2qOPDcUnObxd+eXxa4faAJGliesnRWgA0B4IrAAC4YXqxvNLGimtrevi+HsVCQZ0ady64Wmt1KpnW8eG49u3pZMUVgCQp4HQBAACgeUwurEoiuLYqn89oJBHXqWRapZKVz2caXsN7k4uazub0bCKulVxBySmCKwBWXAEAwAbVGa57Ca4t63girtnldb1zbcGR93+5stp7fDiuxGBEF9PLKhRLjtQCoHkQXAEAwA2phXJwHeCMa8s6NtQnY+RYd+FT42l9+L5u9Uc6lOiPaL1Y0qUZOgsDrY7gCgAAbkgtrikU9CvS0eZ0KXBILNyuj+zr0cvj0w1/7+oYnGcrTaKGByOSxHZhAARXAADwI6kFRuGgvF34+1fmNb+y3tD3rY7BeXa4X5J0OB6WMQRXAARXAACwQYoZrpB0fLhfJSu9cj7T0Pd9eXxa3R0Bfez+PZKkzqBfD0S7NDHFSByg1RFcAQDADVMLaxrs7nS6DDjso/t71NPZdqNRUiNUx+AcS8QV8P/oR9TEQETjrLgCLY/gCgAAJEnFktVUNqfBnnanS4HDAn6fnhnq06lkWtbahrznxjE4GyUGwrqUWVauUGxIHQCaE8EVAABIkmaWciqWrAZ7WHGF9GwirnQ2p/cmFxvyfhvH4GyUGIioULJ6P7PckDoANCeCKwAAkCRNVkbhDDIKByo3aJIaNxZn4xicjRID1c7CnHMFWhnBFQAASCo3ZpKkvTRngqT+7g49tLdbpxpwzvXWMTgbHYqH5PcZJVOccwVaGcEVAABIKo/CkaQBVlxR8exwXG9dnlN2LV/X97l1DM5G7QG/DsS6GIkDtDiCKwAAkFRecW3zG8VCQadLQZM4noirULJ69fxMXd/n1jE4txoejBBcgRZHcAUAAJLKK679kQ75fMbpUtAkPvFgr8Ltgbqec73TGJyNhvojujy7orU8nYWBVkVwBQAAksrBlfOt2KjN79PTR2I6NT5dt7E4dxqDs9HwYETWSuenadAEtCqCKwAAkFTeKjxAcMUtjif6dX1hrW6h8U5jcDZKDIQlie3CQAsjuAIAAFlryyuuNGbCLaqB8uU6dRe+0xicjR6MhRT0+zROcAVaFsEVAABocbWg1XxRg6y44hb79nRqqD9cl3OudxuDs1Gb36dD8ZAmmOUKtCyCKwAAuDHDleCKzRxPxPXG+7NaWS/U9Ll3G4Nzq8RAROPMcgVaFsEVAABocmFVkjTIVmFs4tnhfq0XSzpzobZjce41BmejxEBY1+ZXtZSrbXgG4A4EVwAAoKnKiusAwRWb+OTBXnW2+Wu6XXgrY3A2SgxEJEkTnHMFWhLBFQAAaHKB4Io7aw/49dThWE2D61bG4Gz0o+DKOVegFRFcAQCAphbX1BcOKhjgRwNs7vhwXJdnVvR+Zrkmz9vKGJyN7o92qaONzsJAq+JvJwAAoNTCGo2ZcFfHKyujp8ana/K8rYzB2cjvMzrSH2aWK9CiCK4AAECTC2s0ZsJdPRgL6WBfSC/XYLvwVsfg3CoxECG4Ai2K4AoAADS1yIor7u14Iq7XLs5oLV/c1XO2MwZno8RARFOLOS2s5Hf1/gDch+AKAECLW8sXNbeSZ8UV93Q8EddavqQ33p/d1XO2MwZno+FKg6bkNKuuQKshuAIA0OKqo3AGezodrgTN7olDMQUDvhuNlXZiu2NwNhoaCEsS24WBFkRwBQCgxVVH4bDiinvpDPr1+MGoTiV33qBpu2NwNtq3p1OhoF/JFMEVaDUEVwAAWtyPVlwJrri3Z4f7dSG9rCuzKzt6/XbH4GxkjNHQQERJZrkCLYfgCgBAi7ux4kpwxRbcGIuzw+7C2x2Dc6thOgsDLYngCgBAi0strCnSHlC4PeB0KXCBw/GQ9u3p3FFw3ekYnI2GBsKaWV7XzFJux88A4D4EVwAAWlxqYU0DrLZii4wxenY4ru+cz2i9UNrWa3c6Bmej4cFKZ2G2CwMtheAKAECLSy2uaS/BFdtwPBHX8npRb17e3licnY7B2ShRHYnDdmGgpRBcAQBocamFNQ3QURjb8NSRPrX5zba2C98YgzO0/TE4G/VH2tXdESC4Ai2G4AoAQAsrFEtKL+VYccW2hNsDevTBqE5tY57rjTE4uzjfKpW3Kg8P0qAJ8ILtdCcnuAIA0MIyS+sqliwrrti248NxjaWySlW6Ut/Lbsbg3Ko6Esdau+tnAXDOV76d3PK9BFcAAFpYqjrDleCKbaqunI5ucbvwbsfgbDQ8ENHCal7TWToLA261Xijp2+emtnw/wRUAgBaWWliVxAxXbN/wQEQD3e16OTl9z3trMQZno6GBsCQaNAFu9uqFjLJrhS3fT3AFAKCFVbd5ElyxXcYYHU/EdXoio0Lx7mNxajEGZ6PhSmfh8RTBFXCrb72T2tb8cIIrAAAtLLWYU9DvU7Qr6HQpcKFnh/uVXSvoe1fm73pfLcbgbBQLtysWCmqCWa6AKxWKJf39eym98NDW/zGL4AoAQAtLLayqv7tdPp9xuhS40NNH+uT3mbt2F67VGJxbJQYiGmerMOBKr78/q7mVvF58eO+WX0NwBQCghaUW1xiFgx3r6WzTx+7fc9d5rrUag3OrxEBYE1NZOgsDLvTNs5PqbPPreGLrfy4QXAEAaGGphTVG4WBXnh2O651rC0rfocNvLcfgbJQYjGh5vahr86s1fS7cbzyV1flpVuObVbFk9a2zU3ruaFydQf+WX0dwBQCgRVlrWXHFrh1PlM+onZ7YfNW1lmNwNkpUGjRxzhW3+vKfvaVf/YvvO10G7uCty3PKLOW2tU1YIrgCANCyFlbzWsuXWHHFrnz4vm71hYObbheu9RicjRL9lc7CnHPFBpdnlnUhvax3ry9qfmXd6XKwiW+8M6lgwKfnjm6vyzjBFQCAFpVaLI/C2dvT6XAlcDOfz2hkKK7RZFrF0s3nTWs9Bmejnq42DXS3M8sVNzk5Vp4rbK302sVZh6vBrUolq797N6WRofi2RuFIBFcAAFrW5I0Zru0OVwK3Oz4c19xKXu9cW7jpeq3H4NwqMRAhuOImJ8fTejDWpc42v85cyDhdDm7xg6vzmlxY008+Mrjt1xJcAQBoUVM3gisrrtidY0NxGVMOqlX1GoOzUWIgovPTS7et9KI1rawXdObijH78oQE9eqBXZy7OOF0SbvHNsym1+Y1eeGhg268luAIA0KImF9ZkjNQfYcUVuxMNBfWR/TePxanXGJyNhgciWsuXdGV2pW7vAfc4c2FG64WSnj/ar6cO9yk5tXTHbtdoPGutvnl2Uk8d7lNPZ9u2X09wBQCgRU0trqkv3K62Oq2GobUcT8T1gyvzmlsuN8Sp1xicjYYGwpLEdmFIkk6MTSsU9OvRA7168nBMklh1bSLvXl/UldnVHW0TlgiuAAC0rMmFNQ3SURg18uxwXCUrnT5fPldYrzE4Gw1VRuIQXGGt1cvjaT19pE/tAb8evq9bkfaAzlwguDaLb56dlN9n9BMfIrgCAIBtmFpc0yAzXFEjH92/R3u62nRqPF3XMTgbhdsD2renU0lmuba8ieklXZtf1fOVESsBv0+PH4rSoKlJlLcJp/T4waiioeCOnkFwBQCgRbHiilry+4yODcV1KpnW6Yl03cbg3Gp4kM7CKG8TlnTTf3NPHIrp0syKrs+vOlUWKiaml3QxvawXH9m742fcM7gaY75qjJk2xpzdcC1qjPm2MWai8ntv5boxxvyuMea8MeaHxpiPb3jN5yv3TxhjPr/jigEAwK6trhe1sJpnxRU1dTwRV2Ypp3976kJdx+BsNDQQ1sX0svLFUt3fC83r5Ni0PrS3+6Y/05463CdJbBduAt94Z1LGSJ/+8Pa7CVdtZcX1jyV95pZrvybpJWvtkKSXKp9L0ouShiq/viTp96Vy0JX0G5Iel/SYpN+ohl0AANB4qcXKKBxWXFFDI4lyUDh7bbGuY3A2Gh6IaL1Y0uWZ5bq/F5rTwmpeb16e03NHb96afnQwot6uNho0NYFvnU3p0Qd7d3Xm/Z5/mlhrRyXN3nL5s5K+Vvn4a5J+asP1f2/LXpO0xxizV9KnJX3bWjtrrZ2T9G3dHoYBAECDpG7McCW4onb6Ix368H3dklT3861ViRsNmjjn2qpemcioWLJ67pat6T6f0ROHYjpzYUbWMuvXKRfTSxpLZfXiwzvfJizt/IzrgLV2UpIqv1f/K9kn6cqG+65Wrt3pOgAAcEBqsXzmi+CKWnv+aL/8PqPjicYE1yP9YRkjjac459qqTo5Pa09Xmz72wO0bOp88HNO1+VVdmeWcq1O+eTYlSfrMwzvrJlwVqEUxG5hNrtm7XL/9AcZ8SeVtxnrggQdqVxkAALghtZCTxFZh1N4vPXtYn/7woPob9N9WR5tfD0a7NDFNcG1FpZLVy+PTGhmKy++7PXI8VZnn+p0LGT0QI1s44VtnU/ro/Xt0357OXT1npyuuU5UtwKr8Pl25flXS/Rvu2y/p+l2u38Za+wfW2kettY/G4435lzoAAFrN1OKaIh0Bhdpr/W/YaHVdwYAe3tfT0PdMDETYKtyizl5fUGZp/bbzrVWH42HFI+2cc3XIldkVvXNtQT+5y9VWaefB9W8kVTsDf17SX2+4/guV7sJPSFqobCX+O0mfMsb0VpoyfapyDQAAOGByYZXVVnhGYiCi9zPLyhWKTpeCBjsxNi1jpOOJzUcvGWP05KGYvsM5V0d8q7JNeLfnW6WtjcP5fySdkTRsjLlqjPmipN+R9BPGmAlJP1H5XJK+IemipPOS/p2kL0uStXZW0j+X9N3Kr9+sXAMAAA5ILeY43wrPSAxGVCxZvZ+hs3CrOTme1o/dv0fRUPCO9zx5OKZ0NqcLaVblG+2bZyf1ob3deiDWtetn3XN/kLX2v7vDl17Y5F4r6Zfv8JyvSvrqtqoDAAB1kVpYVaKfIznwhsRAWFK5QdPRwW6Hq0GjZJZy+uHVef3PP564633Vc65nLszoSH+kEaVB5e71b38wr//lU3f//2er6j9cCwAANJVCsaR0Nqe9rLjCIw72heT3GU1wzrWlnBpPy1rpuaObbxOueiDapX17OvWdC5xzbaRvnZ2UJH2mBtuEJYIrAAAtJ72UU8lKAwRXeER7wK+DfSGNT9FZuJWcHJ9WPNKuD+29+yq7MeV5rq9dnFGpxDnXRvnm2ZSG+sM60h+uyfMIrgAAtJjUwpokseIKT0kMhDVBcG0ZhWJJo8m0nhuOy7fJGJxbPXU4prmVvMaY99sQ6WxO3700qxcfqc1qq0RwBQCg5VSD6wBdheEhiYGILs+uaHWdzsKt4O0P5rW4VtBzw3ffJlz1ZPWcK2NxGuLv30upZKUXazAGp4rgCgBAi0ktVldcdzcMHmgmiYGIrBWdY1vEibFpBXxGzwz1ben++/Z06kCsS2cuZOpcGaTyGJwDsS4dHaxdMyyCKwAALSa1sKZgwKferjanSwFqJjFQ/gF5nK2gLeHl8Wl98kBUkY6t/zn25OGYXr84q0KxVMfKMLe8ru9cmNGLj+yVMffexr1VBFcAAFpManFNg90dNf2BAnDagViXgn6fktMEV6+7Pr+qsVRWzx3d3kivJw/3KZsr6N3ri3WqDJL07XNTKpZsTbcJSwRXAABazuRCObgCXhLw+3QoHlKSFVfPOzk+LUl6/h5jcG71xKGopMaec10vlPTlP3tLr7XQ2dpvnU1p355OPbKvp6bPJbgCANBiphbXNEhHYXhQYiCiJLNcPe/kWFr7ezt1OL69MSv9kQ4N9YcbOs/1W++m9I13Uvrtb5yTtd4fxbO4ltcrExm9+PBgzXf1EFwBAGgh1tryiivBFR40PBjRtflVLeUKTpeCOlnLF/Xq+YyeP9q/o2D01OGY3rw0q/VCY865/umZy/L7jH5wdaElOhqfODet9WJJLz5S223CEsEVAICWMr+S13qhxCgceNJQf3kFjnmu3vXG+7NazRe3PAbnVk8ejmllvagfXp2vcWW3G0st6o1Ls/rVF4bUF27Xvz11se7v6bRvnp3UQHe7PnZ/b82fTXAFAKCFTC5UR+EQXOE9w5XRG0mCq2edHJ9We8CnJw7FdvT6xw/GZIwasl34T1+7rGDAp5974kF94ekDGk2m9e71hbq/r1OWcwW9PJ7WZz48KJ+v9s3/CK4AALSQqcoMV1Zc4UX393apo83HOVcPOzk2racOx9QZ9O/o9b2hoB4a7NaZOgfX7Fpe/+nta/qHH7lPvaGgfu6JBxVuD3h61fXl8bRyhZI+8/Deujyf4AoAQAtJLbLiCu/y+YyG+iOsuHrU+5llXZpZ0XPb7CZ8q6cOx/TWB3NayxdrVNnt/vP3rml5vaiff/JBSVJPZ5v++8cf0P/3w+v6YGalbu/rpG+enVQsFNRjB6N1eT7BFQCAFjK5sCZjpHik3elSgLoYGggTXD3qxFh5DM5Oz7dWPXk4pvVCSW9/MFeLsm5jrdWfvHZZj+zr0Uf3/2gkzC8+fVB+n9G/O+29Vde1fFEnx6b1qQ8Pyl+HbcISwRUAgJYytbCmeLhdbX5+BIA3DQ9ENLWY08JK3ulSUGMvj0/rSH9Y90e7dvWcxw5G5feZum0XfuP9WSWnlvTzTzx4U+fjwZ4O/dcf26evv3lFmaVcXd7bKaPJtJbXi3rx4dp3E67iby0AAFrIJDNc4XGJgUqDpmnvrrquF0r6xjuTDRvp0gyWcwW9fnFWzw3Hd/2sSEebHt7XU7fg+ievXVZ3R0D/8KP33fa1L438/+3deXjU5b338c+dPSEhIXvCmoSwCQEksoqAWkUFt6q1x7V1qdX21PpoH2vPqa2P9jndbXtOca+2WtfWBVxRQVAWAdkCYScJkExCgEzWyTb3+WMSGhUwhNkyeb+uywtmMvx+38gNzGfu5Zunlna3nllR4pN7B8o7RQ4lxkZqWl7PDs3qDoIrAAB9SKXTpUwOZkIIG9FxsvB2R2gG19Z2t77//Ge6/bnP9Nzq0kCX4zef7KpWS7v7lPe3dpqel6IN+2rU4OWev1V1Lr1T5NCVhYOPeYDU8PR4nTcmQ8+sKAmZfsMtbW4tLq7U18Zk+HQ1D8EVAIA+pMLZxIwrQlp2YozioyNCspdra7tb//78er27pVIJMRFatKki0CX5zZLtBxUfHaHCod45+Gdabora3FZrS727z/XFT/epzW11zZQhx33NbbPyVOtq0wuflnn13oHyye5q1bnafLpMWCK4AgDQZzS2tKnW1UZwRUgzxig/I17bQyy4trW7decLG/R2kUM/nTdGt83K07rSIyqvaQp0aT5nrdXS7VWamZ+qqAjvxJfCYQMUGW60Yne1V64neX6P/v5pmWbmpyo3Lf64r5s4ZICm5CTrieV7Q2K59zubHYqPjtCZ+ak+vQ/BFQCAPsLh9LTCYakwQt2I9ATtDKFerm3tbt354ga9ublC/3HRaH37zBzNK/D0ynyzD8y6bnPUqcLpOuXThLuKi4rQxMEDtMqL+1w/2FalCqdL104d+pWvvW12nhy1Lr2+4YDX7h8Ibe1uvbfVoXNGpys6ome9dbuL4AoAQB/R2cOVGVeEuhGZCTrU0BISJ7e2tbt110sbtWhThe67cJRunpkrSRqa0k8FgxK1cFN5gCv0vSXbPW1wZnvhYKaupualaPMBp5xN3jmB+tlVpcpKjNE53diHO3tEmkZlJujRZXvkdluv3D8QVu89rCONrT5fJiwRXAEA6DOYcUVfMSLDs0yzt/dzbXdb3f3yRr2xsVz3XjBKt56V97mvzyvI0qb9TpUeaghQhf6xZFuVxg7sr3Qv/901PZ4K0XEAABucSURBVC9FbutpX3Oq9hys1/Kd1fq3yUMU0Y0Diowx+u7sPO2qqtf7xZWnfP9AebuoQrGR4Zo1wnuz4cdDcAUAoI9gxhV9xciOlji9eblwu9vqnpc36rUN5brn/JG6bVbel15zUYGn3UooH9LkbGzVutIjXl0m3GnikCRFR4R5pS3Oc6vLFBFm9I3Jg7v9ay4al6VBA2L1yEe7ZW3vm3Vtd1u9u6VSc0alHfMEZW8juAIA0Ec4nC71j4lQXFREoEsBfCotIVqJsZG99oCmdrfVj17ZpH+uP6C7zxuhO+YMP+brBibFatLQAVq4MXSXC3+086DcVl5rg9NVdES4CocNOOUDmppa2vXy2n2aOzZT6Qnd/2AwIjxMt8zM1WdlNVpT4t3Tjf1hXekRHaxr1tyxWX65H8EVAIA+wuF0KSsxNtBlAD5njNHIjIRe2RLH7ba69x+b9I/P9uuH547Q987OP+Hr5xdkaZujTruqet/32h1Lt1VpQFykxg9K8sn1p+WmaJujTocbWnp8jYUby1XratN13TiU6YuuKhys5H5ReuSj3T2+f6C8XVShqIgwne2DDxWOheAKAEAf4ah1KYNlwugj8jPitd1R16uWYLrdVve9ulkvr9uvH5yTrx+ce+LQKkkXjsuSMdLCjaG3XNjttlq646BmjUhTeJjxyT2m5XlauKza07PlwtZa/XVViUZkxGtyzsn3mI2NCteN04fpw21V2uao7VENgeB2W71T5NBZ+WmKj/bPKh6CKwAAfYTD6VJm/+hAlwH4xcjMBNW62lRV1ztOFna7rX7yWpFeWLNP3z97uO7sRmiVpPT+MZqSk6xFm8p7VUjvjo37a3S4ocUny4Q7FQxKVFxUeI/3uW7c71TRgVpdN3WojOlZuL5+2lDFRYXr0Y/29OjXB8LG/TWqcLr8cppwJ4IrAAB9QGu7Wwfrm5XJUmH0EfnpngOatjuCfwmttVb/+XqRnv+0THfMydNdXxtxUiFo/vhs7T7YoOKK4P9eT8aS7QcVZqRZI7zbBqeryPAwTc5J7vE+17+tLFW/qHBdOnFgj2tIiovS1WcM0Rsby7X/SGOPr+Mvbe1u/fa9HYqKCNO5ozP8dl+CKwAAfcDBumZZSysc9B29pSWOtVY/fX2Lnltdpttm5enu80ae9MzdBWOzFB5mtCjEerou2Val04cMUFJclE/vMz0vRbsPNqiq4+T17jrS0KKFm8p12ekDlRATeUo13DwzR0bSE8v3ntJ1/OGht4r18a5qPXjJWCXGndr3fTIIrgAA9AEVHT1cs9jjij4iJT5aqfFRQR1crbX6+cKt+tuqUn3nrFz937knH1olKblflGYMT9XCEFouXFXn0uYDTp8uE+40Ldezz3XlSe5zfXndPrW0uXVtDw5l+qLspFhdPCFbL67ZpyOncFCUr724pkx/+aRE356Ro6vO6H7rH28guAIA0AdUdswkZDDjij4kPz1BO4K0l6u1Vg8s2qqnV5To5jNzdO8Fo3q8R1KS5hVkad/hJm3a7/RilYGzdPtBSdLskb5bJtxpTHZ/9Y+JOKl9rm631bOryjR5WLJGZfb3Sh23zcpTU2u7nllZ4pXredvaksP6j9eKNDM/VfddOMrv9ye4AgDQBziYcUUfNDLT0xIn2GYhrbV68M3iozNXP7lo9CmFVkk6/7RMRYabkOnpunR7lTL6R2tMlndC4YmEhxlNyU3RipMIrst2HlTZ4UZdO+3UZ1s7jchI0Lmj0/X0ihI1trR57brecKCmSbc9u06DBsTpv795uiLC/R8jCa4AAPQBjlqXoiLClOTH/UhAoOVnxKuhpV0HapoCXcpR1lr9/7e36cmP9+rG6cP0n/NOPbRKUmJspGaNSNObmyvkdgdXUD9Zre1uLd9RrTkj073y/6Y7puelqOxwY7cPR3p2ValS46M09zTvnqp726w81TS26sU1+7x63VPR2NKmW55Zq+ZWtx6/vtCv+1q7IrgCANAHOJwuZSXG+O1NIBAMRmZ4ThYOln2u7W6r/3p7mx5btkfXTxuq++eP8eqfyXkF2apwurSu7IjXrhkIa0uOqK65zS/7WztNy0uRpG4tF953uFEfbKvS1WcMUVSEd+NU4bBkFQ4doCeW71Vru9ur1+4Ja63ueXmTih21+uO/TdTw9PiA1UJwBQCgD3A4XexvRZ+TfzS4Bn6f697qBn3j0ZV6dNkeXTt1iH5+8Wle/yDp3DEZio4I06Jevlx46fYqRYYbzRie6rd7jkhPUEq/qG4F1+c/LZOR9M0pQ3xSy22z8nSgpikoTon+04e79ObmCv34glGaM9J/HyQcC8EVAIA+wFHrYn8r+pzE2Ehl9o/R+rIjAdvn2u62emL5Hs19eJl2VNbpt1eO1/+7ZKxPVj/ER0fo7FHpenOzQ+29eLnwh9uqNDknWfHREX67Z1iY0dTcFK3cc+iEY6W5rV0vrtmnc0ZnaGCSb/pinz0qXSMy4vXI0j0B3Z/9TpFDv1u8Q5dPHKhbZuYGrI5OBFcAAEKctVaOWhc9XNEnzR6Zpne3VOqyP6/QutLDfr13SXWDrn5spR58s1gzhqdq8V2z9PVJg3y6ZH/++GxV1zdr9Um2dgkW+w43amdVfUBm96blpajC6VLJoePvc32nyKFDDS26zgstcI4nLMzoO2flaXtlnZZsr/LZfU6kuKJWd720QeMHJ+kXl48Lim0mBFcAAELckcZWtbS5lcmMK/qghy4bp19fUaDymiZ9fcFK3fHcZyo7QTDxBrfb6qmP92ruH5Zpu8Mzy/rkDYV+Wa4/Z2S64qLCtXBThc/v5QtLd3ja4Phzf2un6d3Y5/q3laUalhKnM328jPniCdnKTozRI0v3+PQ+x3Kovlk3P7NWCTERevy6SYqJDPd7DcdCcAUAIMRVOD0nqjLjir4oPMzoysLBWnrPbN15br4+3Falc3/3kX7xVrGcTa1ev59nlnWVHli0VdPzUvXeD30/y9pVbFS4vjYmQ28XVQTF4T4na8m2Kg1NiVNuaj+/3zsntZ8y+kdrxe7qY359a3mt1pYe0bVThyoszLe/n5HhYbppZq4+LTmsdaX+O2yrpc2t7z73marrm/XYdYVKD6J/NwiuAACEuMpaTw9XZlzRl8VFRejOc0doyd2zdfGEbD2+fI9m/3qJ/rqyxCsBr+ssa7GjVr++okBP3lAYkD938wqyVdPYqk92HTuABStXa7tW7PZvG5yujDGanpeqVcfZ5/rs6lJFR4TpikmD/FLP1WcMVmJspB75aLdf7met1c8WbtGnew/rV1cUaPzgJL/ct7sIrgAAhLgKJ8EV6JSZGKPfXDleC793pkZl9tdPX9+i8x9epg+KK3t8EE7poQZd/bhnlnVabooW/3CWriwcHLB9gWeNSFVCTIQWbuxdy4VX7jkkV6tbs0emBayGabkpqq5v0c6qz59EXetq1WvrD+ji8dlKiovySy39oiN0w/RhWry1UruqfN/S6dlVpfr76jJ9d3aeLpkw0Of3O1kEVwAAQlyl06UwI6XFRwe6FCBojB2YqL/fMkVPXF8oWemmZ9bqmidWa0u5s9vXcLutnv5kr+Y+vFzFFZ5Z1qduPCPgHxJFR4Tr/NMy9d4Wh5rb2gNaS3eVHmrQHz/YqZjIME3NTQlYHcfr5/rqZwfU2NKu66b57lCmY7lx+jDFRIbpkY98u9d1xa5q/WzhVp0zKl13nzfSp/fqKYIrAAAhrsLpUlpCtCLC+Wcf6MoYo3PHZOjdH56ln198moorajXvTx/rR69sPLrE/nhKDzXom4+v0s8WbtWU3GS998OzAjrL+kXzCrJU19ymj7YfDHQpJ9Ta7taCpbt13u+XaWdlvX5x2biAHgY0ODlOgwbEfm6fq7VWf1tVqvGDElUwyL/LZ5P7RekbhYP1+oYDR88r8LayQ426/e+fKSe1nx6+eoLCfbx/t6f4FwwAgBBHKxzgxCLDw3TD9GFaes8c3TIzV6+uP6DZv16qP7y/U40tbZ97bddZ1q3ltfrVFQX6y41nKCvRNz09e2rG8FQNiIvUoiA+XXh92RHN/9PH+uU72zRnZLrev2uWLj/dP/tHT2R6XopW7Tksd0cv3FV7DmtXVb2u9WELnBO5eWau3FZ6cvler1+7vrlNN/91jayVnri+UAkxkV6/h7cQXAEACHEOpyvgSxeB3iAxNlL3XTha7981S3NGpen37+/QnN8s1Svr9svttio71Hh0lnVyTrLeu+ssXRVEs6xdRYaHae7YLL1fXKmmluBaLlznatX9rxfp8gUrVNPYqseum6RHrpsUNH9PTctLkbOpVVsraiV59n4mxkZq/vjsgNQzODlO8wqy9PynZfpox0HVurxzGrbbbXXnCxu0+2CD/nzN6RoWgJOcT0ZEoAsAAAC+5ah1He1PCOCrDU3ppz9fM0lrSg7rwTeLdffLG/X4sj3ad6RR4cboV18v0JWF/mtx01Pzx3vCzofbqnRRQVagy5EkvbvFoftf36LKOpdumDZMd58/UvHRwRVJpuV6erSu3H1IaQnReneLQ9+aMSygS5hvnz1c72+t1A1PfSpjpOFp8Zo4JEkThwzQxCFJyk9POOklvr9bvEPvF1fqZ/PHaIaP+9J6Q3CNEgAA4FUNzW2qc7UpM8iWMQK9wRnDkvXqd6dr4aZy/ea97ZqSk6yHLhun7KTe8edpSk6K0hKitWhTecCDa4WzSfe/vkXvba3U6Kz+euS6SZoQZO1WOmUmxig3tZ9W7jmkxpZ2tbmtrpkSmGXCnUZmJmjlfedo0z6n1pcd0fp9NVq8tVIvrd0vSYqPjtD4wYmaONgTZCcMTlLKCQ7ke2Njuf57yS5dfcZg3TB9mJ++i1NDcAUAIIQ5jvZw5URhoCfCwowumTAwKNuDfJXwMKMLx2bqhTX7VN/cFpCZzXa31bOrSvXrd7erze3WvReM0k1n5igyyA+Lm5aXotc3lGtLuVNnjUgLimW0/WMidWZ+qs7M98yOWmtVcqjRE2TLarR+3xEt+Gi32jv25g5NidPEwf+alR2d1V+R4WHavN+pe17eqDOGDdADl4wN+pUDnQiuAACEsMrOHq79e8cMEQDvmj8+W8+sLNX7Wyt16UT/hu/iilr9+J+btWFfjWbmp+qhS8dpSEqcX2voqel5qXpudZnqm9v04KWBnW09HmOMclL7KSe139FDrZpa2rX5gPNomF2x+5Be21AuSYqOCFPBoESVHmpUany0Flw7SVERwf0BQlcEVwAAQlhFZ3ANkkNPAPjX6UMGKCsxRgs3lvstuDa1tOsPH+zUE8v3KDE2Un+4eoIuHp/da2b2JGlqbrIkaWBSrM4elR7garovNipck3OSNTnHU7+1VhVOl2dGtmOJcWR4mB67fpJSe1lvb4IrAAAh7OhSYdrhAH1SWJjRvIIsPb2iRM7GViXG+bbdyfKdB/WTV4tUdrhRVxUO0n0XjlZSXJRP7+kLKfHRunH6ME0ckhS0fU27wxij7KRYZSfFBnyf86kiuAIAEMIcTpcSYyMVGxW40zABBNa8gmw9vnyv3t3i0FVnDPbJPQ7VN+vBN4v16voDyk3tp+dvmappvfw0859dfFqgS0AXBFcAAEKYo9alLJYJA31awaBEDUmO08JN5T4Jrou3VuqeVzaqoblN/35Ovm6fnRfQ1jEITQRXAABCmMPpUgbLhIE+zRjPcuFHl+3RofrmE7ZJOVl/+WSvHli0VeMGJuq3V45XfkaC164NdNV7jpECAAAnjRlXAJLndOF2t9XbRQ6vXK/dbfXAwq36+cKt+troDL146zRCK3yK4AoAQIhqbXerur6ZGVcAGpWZoLy0flq0qfyUr9XU0q7bn1unpz7Zq2/NGKYF105iHz18juAKAECIqqprlrVixhVAx3LhbK3ee1iVHaeN90R1fbO++fgqvbe1Uj+dN0b3zz+tV5+6i96D4AoAQIhyOJskSRkEVwCS5o/PkrXSW5srevTr9xys1+V/XqHiilotuGaSvn1mjpcrBI6P4AoAQIhyOJslMeMKwGN4eoJGZSZo4caTXy68puSwLl+wQg3NbXrh1qmaOzbTBxUCx0dwBQAgRFV0zLhmsscVQIf547P1WVmNDtQ0dfvXLNxYrmueWK3kuCj98/bpmjhkgA8rBI6N4AoAQIiqrHUpOiJMibGRgS4FQJCYX5AtSXqzG4c0WWu1YOluff/59Ro/KFH/+O50DU3p5+sSgWMiuAIAEKIqnJ5WOMZwcAoAjyEpcRo/KFELN554n2tbu1v/8VqRfvnONs0ryNLfbpqiAf2i/FQl8GUEVwAAQlRlrYtWOAC+ZF5BtjYfcKqkuuGYX29obtMtf12r51aX6bZZefrj1RMVE0m7GwQWwRUAgBDVOeMKAF1dVJAlScfs6VpZ69JVj67URzsO6qHLxureC0YpjHY3CAIEVwAAQpC1VlW1zbTCAfAl2UmxKhw6QIs2fX658HZHnS77n0+0t7pBT95whq6ZMjRAFQJfRnAFACAEHW5oUUu7W1ksFQZwDPMKsrTNUaedlXWSpE92VeuKBSvU6rZ66TvTNGdUeoArBD6P4AoAQAiqcLokSZnMuAI4hgsLshRmpIWbKvTKuv264alPlZUUo9fumKGxAxMDXR7wJRGBLgAAAHhfZW1ncI0NcCUAglF6Qoym5KToyeV71NDSrul5KVpw7STaZyFoMeMKAEAIOjrjylJhAMdx6cRsNbS06/LTB+rpb00mtCKoMeMKAEAIqqx1KTzMKC0hOtClAAhSV04arNFZ/TVuYCL9nhH0CK4AAISgCqdLafHRCqeNBYDjCAszKhiUFOgygG5hqTAAACGostbFwUwAgJBBcAUAIARVOF3sbwUAhAyCKwAAIajSyYwrACB0EFwBAAgx9c1tqmtuI7gCAEIGwRUAgBDj6GiFk0VwBQCECIIrAAAhpjO4ZrDHFQAQIgiuAACEGEetJ7hyOBMAIFQQXAEACDEOZ5MksccVABAyCK4AAIQYR61LSXGRiokMD3QpAAB4BcEVAIAQ46CHKwAgxBBcAQAIMY5aergCAEILwRUAgCBhrfXKdRxOF61wAAAhJSLQBQAA0Jc1trRp6faDervIoSXbqtTmdispNkpJcZFKjI3s8mPU0cdJsVFf+Fqk4qMjZIxRS5tb1fUttMIBAIQUgisAAH5W52rVh9uq9PZmh5buqJKr1a3kflG6cFymEmMj5WxqVU1jq2qaWlVS3aiaphY5m1rlanUf95rhYUaJsZ4AK4kZVwBASCG4AgDgB87GVi0urtQ7RRVatqNaLe1upSVE68pJg3XBuExNHpasiPAT7+Bxtbb/K9Q2esJsTVOrnI2tqmlqUU1jq5xNrRqZmaAZw1P99J0BAOB7fg+uxpi5kv4gKVzSE9ba//J3DQAA+MOh+mYt3lqpt4ocWrGrWm1uq+zEGF07daguHJep04cMUFiY6fb1YiLDFRMZzjJgAECf49fgaowJl/Q/kr4mab+kNcaYN6y1W/1ZBwAAvlJV69K7Wxx6u8ihVXsOyW2lIclxumlmji4Ym6XxgxJlTPfDKgAA8P+M62RJu6y1eyTJGPOCpEskHTO4NrW2a0u504/lAQBw8txuaU3JYb1dVKG1pUdkrZSX1k93zBmuuWMzNSarP2EVAIBT4O/gOlDSvi6P90uacrwX76qq10V//NjnRQEA4A2jMhN05zkjdOG4TOVnJAS6HAAAQoa/g+uxPm7+XNM6Y8ytkm6VpIxBw/TodZP8URcAAKdkREaCclL7BboMAABCkr+D635Jg7s8HiSpvOsLrLWPSXpMkgoLC+35p2X6rzoAAAAAQNA58bn73rdGUr4xJscYEyXpaklv+LkGAAAAAEAv4tcZV2ttmzHme5LelacdzlPW2i3+rAEAAAAA0Lv4vY+rtfYtSW/5+74AAAAAgN7J30uFAQAAAAA4KQRXAAAAAEBQI7gCAAAAAIIawRUAAAAAENQIrgAAAACAoEZwBQAAAAAENYIrAAAAACCoEVwBAAAAAEGN4AoAAAAACGoEVwAAAABAUCO4AgAAAACCGsEVAAAAABDUCK4AAAAAgKBGcAUAAAAABDWCKwAAAAAgqBFcAQAAAABBjeAKAAAAAAhqxlob6BqOyxhzUFJpoOtAn5YqqTrQRQBewFhGKGAcI1QwlhEKvDWOh1pr077qRUEdXIFAM8astdYWBroO4FQxlhEKGMcIFYxlhAJ/j2OWCgMAAAAAghrBFQAAAAAQ1AiuwIk9FugCAC9hLCMUMI4RKhjLCAV+HcfscQUAAAAABDVmXAEAAAAAQY3gCnQwxjxljKkyxhR1eS7ZGLPYGLOz48cBgawR+CrGmMHGmCXGmGJjzBZjzA86nmcso1cxxsQYYz41xmzsGMs/73g+xxizumMsv2iMiQp0rcBXMcaEG2PWG2MWdTxmHKPXMcaUGGM2G2M2GGPWdjznt/cXBFfgX56WNPcLz90r6QNrbb6kDzoeA8GsTdL/sdaOljRV0h3GmDFiLKP3aZZ0trV2vKQJkuYaY6ZK+qWk33eM5SOSbgpgjUB3/UBScZfHjGP0VnOstRO6tMHx2/sLgivQwVq7TNLhLzx9iaRnOn7+jKRL/VoUcJKstRXW2s86fl4nzxulgWIso5exHvUdDyM7/rOSzpb0SsfzjGUEPWPMIEkXSXqi47ER4xihw2/vLwiuwIllWGsrJE8gkJQe4HqAbjPGDJM0UdJqMZbRC3Usr9wgqUrSYkm7JdVYa9s6XrJfng9mgGD2sKQfSXJ3PE4R4xi9k5X0njFmnTHm1o7n/Pb+IsJXFwYABI4xJl7SPyTdaa2t9XzAD/Qu1tp2SROMMUmSXpU0+lgv829VQPcZY+ZJqrLWrjPGzO58+hgvZRyjN5hhrS03xqRLWmyM2ebPmzPjCpxYpTEmS5I6fqwKcD3AVzLGRMoTWp+z1v6z42nGMnota22NpKXy7NtOMsZ0fvA+SFJ5oOoCumGGpIuNMSWSXpBnifDDYhyjF7LWlnf8WCXPh4mT5cf3FwRX4MTekHRDx89vkPR6AGsBvlLH3qknJRVba3/X5UuMZfQqxpi0jplWGWNiJZ0rz57tJZKu6HgZYxlBzVr7Y2vtIGvtMElXS/rQWnuNGMfoZYwx/YwxCZ0/l3SepCL58f2FsZaVCYAkGWOelzRbUqqkSkn3S3pN0kuShkgqk3SltfaLBzgBQcMYc6ak5ZI261/7qe6TZ58rYxm9hjGmQJ6DPsLl+aD9JWvtA8aYXHlmrpIlrZd0rbW2OXCVAt3TsVT4bmvtPMYxepuOMftqx8MISX+31j5kjEmRn95fEFwBAAAAAEGNpcIAAAAAgKBGcAUAAAAABDWCKwAAAAAgqBFcAQAAAABBjeAKAAAAAAhqBFcAAAAAQFAjuAIAAAAAghrBFQAAAAAQ1P4XzvIWipbHREIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x648 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='UTF-8'), 'fastq')\n",
    "n_cnt = defaultdict(int)\n",
    "for rec in recs:\n",
    "    for i, letter in enumerate(rec.seq):\n",
    "        pos = i + 1\n",
    "        if letter == 'N':\n",
    "            n_cnt[pos] += 1\n",
    "seq_len = max(n_cnt.keys())\n",
    "positions = range(1, seq_len + 1)\n",
    "fig, ax = plt.subplots(figsize=(16,9))\n",
    "ax.plot(positions, [n_cnt[x] for x in positions])\n",
    "ax.set_xlim(1, seq_len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: 0.40 52229\n",
      "1: 1.52 200558\n",
      "2: 3.77 498679\n",
      "3: 4.04 533458\n",
      "4: 4.77 630923\n",
      "5: 4.88 645266\n",
      "6: 2.50 330834\n",
      "7: 2.51 331743\n",
      "8: 2.53 334410\n",
      "9: 2.51 332259\n",
      "10: 4.95 654154\n",
      "11: 2.41 318303\n",
      "12: 2.35 309918\n",
      "13: 2.28 301033\n",
      "14: 2.20 291341\n",
      "15: 2.12 280719\n",
      "16: 2.05 270431\n",
      "17: 1.97 259779\n",
      "18: 1.88 248982\n",
      "19: 1.81 239621\n",
      "20: 1.73 228923\n",
      "21: 1.66 219602\n",
      "22: 1.59 209905\n",
      "23: 1.52 201164\n",
      "24: 1.46 193259\n",
      "25: 1.40 184846\n",
      "26: 1.33 176263\n",
      "27: 1.28 168902\n",
      "28: 1.23 162226\n",
      "29: 1.17 154892\n",
      "30: 1.13 149449\n",
      "31: 1.08 142464\n",
      "32: 1.03 136763\n",
      "33: 0.99 131291\n",
      "34: 0.95 125624\n",
      "35: 0.91 120704\n",
      "36: 0.88 115701\n",
      "37: 0.84 111179\n",
      "38: 0.80 106290\n",
      "39: 0.78 102568\n",
      "40: 22.76 3007221\n"
     ]
    }
   ],
   "source": [
    "recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')\n",
    "cnt_qual = defaultdict(int)\n",
    "for rec in recs:\n",
    "    for i, qual in enumerate(rec.letter_annotations['phred_quality']):\n",
    "        if i < 25:\n",
    "            continue\n",
    "        cnt_qual[qual] += 1\n",
    "tot = sum(cnt_qual.values())\n",
    "for qual, cnt in cnt_qual.items():\n",
    "    print('%d: %.2f %d' % (qual, 100. * cnt / tot, cnt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAIMCAYAAADvmRGtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+UpWdBJ/jvk3Q3kF+dpLvTSSU0cSsgup41zvaAmoODqCEwQtDFHTjqcEaYjKXsqGM5CMPsqASQsRH3zNFiIyDZFUUWkE4QgYhkQWQbG40h0A3kKglQ+dHdSf9KAv3r2T/ubexJ1a2q2133uberP59z+lTd+77Pfb7nVtXT9a33fe8ttdYAAABAK2eNOgAAAABnFkUUAACAphRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmVrWcbP369fXKK69sOSUAAACNfPazn91da92w2H5Ni+iVV16Z7du3t5wSAACARkop9yxlP6fmAgAA0JQiCgAAQFOKKAAAAE0pogAAADSliAIAANDUkotoKeXsUsrflVI+2Lv9baWUbaWUL5dS/qSUsmZ4MQEAAFgpBjki+gtJdpxw+01J3lJrfWqSh5O8fDmDAQAAsDItqYiWUq5I8i+TvK13uyR5TpL39na5OcmLhhEQAACAlWWpR0R/J8l/THKsd3tdkr211iO9219Lcvl8A0spN5RStpdStu/ateuUwgIAAHD6W7SIllJ+NMmDtdbPnnj3PLvW+cbXWm+qtW6utW7esGHDScYEAABgpVi1hH2uSfLCUsrzkzwxyQXpHiG9sJSyqndU9Ioks8OLCQAAwEqx6BHRWuura61X1FqvTPKSJH9Za/3JJB9P8uLebi9LsnVoKQEAAFgxTuV9RF+V5D+UUu5O95rRty9PJAAAAFaypZya+y211tuT3N77/B+SPGP5IwEAALCSncoRUQAAABiYIgoAAEBTA52aCwAAnF6uvfbab33+0Y9+dIRJxieLHHO1zuKIKAAAAE0pogAAsEKdeJRrvtstjUsWOeYaRZaxOzV3ZmYmnU5n3m2zs7NJkomJiXm3T05OZmpqauhZxiVH6yzjkmOhLOOSo3WWccmxUJZxydE6y7jkOJUs45JjubOMS45TyTIuOZY7y7jkOJUs45JjubOMS45TyTIuOYaRBU4nY1dEF/LYY4+NOkKS8cmRjE8WOeYalyxyzDUuWcYlRzI+WeSYa1yyyDHXuGSRY65xyTIuOWAclVprs8k2b95ct2/fftLjp6enkyRbtmxZrkindY5kfLLIMde4ZJFjrnHJMi45kvHJIsdc45JFjrnGJYscc41LlnHIMd8plqN6UZxxySLHXMuZpZTy2Vrr5sX2c40oAAAATSmiAACwQj3+qNYo3yJkXLLIMdcosiiiAAAANHVavVgRAAAwmFEeaXu8cckix1ytszgiCgAAQFOOiAIAwAp24iuijvoI3LhkkWOu1lkcEQUAAKApRRQAAFaox78/5HzvF9nKuGSRY65RZFFEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAYIV6/NtwjPItQsYlixxzjSKLIgoAAEBTq0YdAAAAGJ5RHml7vHHJIsdcrbM4IgoAAEBTiigAAABNOTUXAABo4tprr/3W56M8LVWOuVpncUQUAACAphRRAABg6E484jbfbTlGk2O+uVtkcWouAACc5mZmZtLpdObdNjs7mySZmJiYd/vk5GSmpqaGlg3mo4gCAMAK9thjj406AsyhiAIAwGluoSOa09PTSZItW7a0igOLco0oAAAATSmiAADA0D3+LUFG9XYlcsw1iiyKKAAAAE25RhQAAGhilEf9TiTHXK2zOCIKAABAU46IAgAATVx77bXf+nycjgaOyjg9H62zOCIKAABAU4ooAAAwdCcecZvv9plmnJ6PUWRxai4AALDizMzMpNPpzLttdnY2STIxMTHv9snJyUxNTQ0tG4ooAABwhnnsscdGHeGMp4gCAAArzkJHNKenp5MkW7ZsaRWHx3GNKAAAAE0pogAAwNA9/i1BRv12JaM2Ts/HKLIoogAAADTlGlEAAKCJM/0o6OON0/PROosjogAAADTliCgAAMAIXHvttd/6fNRHR1tnWfSIaCnliaWUz5RS/r6U8vlSyq/37n9nKeUfSyl39P5dPfS0AAAAnPaWcmruN5M8p9b63UmuTnJdKeV7e9t+pdZ6de/fHUNLCQAAsIKceARyvtstjSLLoqfm1lprkoO9m6t7/+owQwEAAKenmZmZdDqdebfNzs4mSSYmJuZsm5yczNTU1FCzMT6W9GJFpZSzSyl3JHkwyW211m29Ta8vpdxZSnlLKeUJfcbeUErZXkrZvmvXrmWKDQAAnG4ee+yxPPbYY6OOwRhY0osV1VqPJrm6lHJhkj8tpXxXklcnuT/JmiQ3JXlVkt+YZ+xNve3ZvHmzI6kAALCCLXRUc3p6OkmyZcuWVnEYUwO9fUutdW+S25NcV2u9r3Z9M8kfJHnGEPIBAACwwizlVXM39I6EppTypCQ/nGRnKeWy3n0lyYuS3DXMoAAAACvF498iZZRv3zKKLEs5NfeyJDeXUs5Ot7i+p9b6wVLKX5ZSNiQpSe5I8rNDzAkAAMAKsZRXzb0zyffMc/9zhpIIAADgDDDKo6CP1zrLQNeIAgAAwKlSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaWsrbtwAAAHASZmZm0ul05t02OzubJJmYmJh3++TkZKampkaeZblzJIooAADASDz22GOjjvAtrbMoogAAAEOy0JHE6enpJMmWLVvOuCyuEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgKUUUAACAphRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgKUUUAACAphRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgKUUUAACAphRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoatEiWkp5YinlM6WUvy+lfL6U8uu9+7+tlLKtlPLlUsqflFLWDD8uAAAAp7ulHBH9ZpLn1Fq/O8nVSa4rpXxvkjcleUut9alJHk7y8uHFBAAAYKVYtIjWroO9m6t7/2qS5yR5b+/+m5O8aCgJAQAAWFGWdI1oKeXsUsodSR5McluSTpK9tdYjvV2+luTyPmNvKKVsL6Vs37Vr13JkBgAA4DS2pCJaaz1aa706yRVJnpHkO+bbrc/Ym2qtm2utmzds2HDySQEAAFgRBnrV3Frr3iS3J/neJBeWUlb1Nl2RZHZ5owEAALASLeVVczeUUi7sff6kJD+cZEeSjyd5cW+3lyXZOqyQAAAArByrFt8llyW5uZRydrrF9T211g+WUr6Q5N2llBuT/F2Stw8xJwAAACvEokW01npnku+Z5/5/SPd6UQAAAFiyga4RBQAAgFOliAIAANCUIgoAAEBTiigAAABNKaIAAAA0pYgCAADQlCIKAABAU4ooAAAATSmiAAAANKWIAgAA0JQiCgAAQFOKKAAAAE0pogAAADSliAIAANCUIgoAAEBTiigAAABNKaIAAAA0pYgCAADQlCIKAABAU4ooAAAATSmiAAAANKWIAgAA0JQiCgAAQFOKKAAAAE0pogAAADSliAIAANCUIgoAAEBTiigAAABNrRrFpDMzM+l0OgOPOz5menp64LGTk5OZmppalizjkmMYWcYlBwAAsHKNpIh2Op3cvWNHnrL24oHGrTlakySHZx8YaNw9+x5aJMsXsmnt+QPkOJIkOTT71YFy3LvvwII5vvyFz2XT2tUDPebxLN/8+s4Bsxzum+NLO+7MxNoy0OOd3fvaHJz93EDjZvfVgfYHAABOfyMpoknylLUX57XPem6TuW785EcW3L5p7fl57Q88c/g5PrFtkRyr86prLhl6jiR506ce7LttYm3Jzz5rTZMcb/3koSbzAAAA48M1ogAAADSliAIAANCUIgoAAEBTiigAAABNKaIAAAA0pYgCAADQlCIKAABAU4ooAAAATa0adQDoZ2ZmJp1OZ6Axx/efnp4eeL7JyclMTU0tS45hZBmXHAAAcKoUUcZWp9PJF3fcmY0XliWPOetYTZLsve9zA831wN66YI6dO+7M+gsHesjkWPfD7vvuHGjY7r39c+zYcWcuumjAGL0c998/WI6HHx5sHgAAWCpFlLG28cKSn/zB4X+bvuvjRxbcvv7C5PofOnvoOZJk68eO9t120UXJD/9Ikxj5i9vazAMAwJnHNaIAAAA0pYgCAADQlCIKAABAU4ooAAAATSmiAAAANKWIAgAA0NSiRbSU8uRSysdLKTtKKZ8vpfxC7/5fK6V8vZRyR+/f84cfFwAAgNPdUt6g8UiSX661/m0p5fwkny2lHH+HwbfUWrcMLx4AAAArzaJFtNZ6X5L7ep8fKKXsSHL5sIMB42tmZiadTmegMcf3n56eHni+ycnJTE1NLUuOYWQZlxwAAKeLpRwR/ZZSypVJvifJtiTXJHllKeVfJ9me7lHTh+cZc0OSG5Jk06ZNpxgXGAedTidf2HFnLrh46WOO1O7Hrz1w50Bz7X9o4Ryf33lnzlk30EPmUC/LP+4aLMuje/rnuHPnncn6wXIcd+fuwXJk98nNAwAwLpZcREsp5yV5X5JfrLXuL6XMJHldktr7+OYkP/P4cbXWm5LclCSbN2+uyxEaGL0LLk6+73ll6PN8+s8XXjbOWZc8/UfbvO7azg8e679xfXLWjw/0t72Tduz9R5rMAwAwLEv67a2UsjrdEvquWuv7k6TW+kCt9Wit9ViS30/yjOHFBAAAYKVYyqvmliRvT7Kj1vrbJ9x/2Qm7/ViSu5Y/HgAAACvNUs4juybJTyf5XCnljt59r0ny0lLK1ememvuVJP9uKAkBAABYUZbyqrl/lWS+C8E+tPxxAAAAWOnavMIHAAAA9CiiAAAANKWIAgAA0JQiCgAAQFOKKAAAAE0t5e1bADhNzMzMpNPpDDzu+Jjp6emBxk1OTmZqamrkORbKAgCMH0UUYAXpdDq5c+fnk/VPGnDkoSTJnbv/YelDdj+2SI4vpKw7f6AUtR5Jknxu11cHG7fnwED7AwCjpYgCrDTrn5RVL3r60Kc58oGdC24v687PquufOfQcSXJk67Ym8wAAy8M1ogAAADSliAIAANCUIgoAAEBTiigAAABNKaIAAAA0pYgCAADQlCIKAABAU4ooAAAATa0adQAAGJaZmZl0Op2Bxx0fMz09PfDYycnJTE1NjTxLvxwAMA4UUQBWrE6nkzt37khZd+FA42o9liT53K77Bhu3Z+8Sslw8YJbay/LAADkeGmgOAGhNEQVgRSvrLsyqFz67yVxHbrl9kSwXZ/ULrht6jsO3fnjocwDAqXCNKAAAAE0pogAAADSliAIAANCUIgoAAEBTiigAAABNKaIAAAA0pYgCAADQlCIKAABAU4ooAAAATa0adQAAoJ2ZmZl0Op2Bxx0fMz09PfDYycnJTE1NjWUOAEZDEQWAM0in08mdO3emrFs/0Lhaux8/t2v3YOP2zL9/N8cXc9a6SwZ6vGO1JEnu2vXwYOP2PDjQ/gAMlyIKAGeYsm591rzg+iZzHbp1a99tZ627JE98wUua5PjGre9uMg8AS+MaUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgKUUUAACAplaNOgAAwCjNzMyk0+kMPO74mOnp6YHGTU5OZmpqauQ5FsoCMGyKKABwRut0Orlz55dy9rpLBxp3rHZPLPv8rv1LHnN0z/0L5vjczi/n7HVPHijH0bo6SfKFXd8YbNyerw60P8ByUkQBgDPe2esuzTkvfMXQ53n0lrctkuPJueD6Xx56jiTZv/XNTeYBmI9rRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgqUWLaCnlyaWUj5dSdpRSPl9K+YXe/ReXUm4rpXy59/Gi4ccFAADgdLeUI6JHkvxyrfU7knxvkp8vpXxnkl9N8rFa61OTfKx3GwAAABa0arEdaq33Jbmv9/mBUsqOJJcnuT7Js3u73Zzk9iSvGkpKAACamJmZSafTGXjc8THT09MDj52cnMzU1NTA44DT16JF9ESllCuTfE+SbUk29kpqaq33lVIu6TPmhiQ3JMmmTZtOJSsAAEPW6XTyhZ135/z1g/3ediRrkiRf3X1ooHEHdt870P7AyrDkIlpKOS/J+5L8Yq11fyllSeNqrTcluSlJNm/eXE8mJAAA7Zy/flOeef1rmsy1besbmswDjJclvWpuKWV1uiX0XbXW9/fufqCUcllv+2VJHhxORAAAAFaSpbxqbkny9iQ7aq2/fcKmW5K8rPf5y5JsXf54AAAArDRLOTX3miQ/neRzpZQ7eve9JslvJnlPKeXlSe5N8hPDiQgAAMBKspRXzf2rJP0uCP2h5Y0DAADASreka0QBAABguSiiAAAANKWIAgAA0JQiCgAAQFOKKAAAAE0t5e1bAACguZmZmXQ6nYHHHR8zPT090LjJyclMTU2NPMdCWWClUEQBABhLnU4nO3fenXXrnjLQuFrXJEl27Tq85DF79tyzYI4v7rg7Gy8eLMdZvRx7H1h6jiR54KH+WWClUEQBABhb69Y9JS984X8e+jy33PK6BbdvvPgp+annvnboOZLkDz9yY5N5YJRcIwoAAEBTiigAAABNKaIAAAA0pYgCAADQlCIKAABAU4ooAAAATSmiAAAANKWIAgAA0NSqUQcAAAAWNzMzk06nM/C442Omp6cHHjs5OZmpqamBx8FiFFEAADgNdDqdfGnH3bl87aaBxq06uiZJ8sjsoYHGfX3fvQPtD4NQRAEA4DRx+dpN+blnv7bJXL93+41N5uHM5BpRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaWjWKSWdnZ/Po3n258ZMfaTLfPXsfyjk52jfLI3sP5MZPbGuQ40DOzWz/HPsO502fenDoOZLk3n2Hc26ZPwsAACxkZmYmnU5n4HHHx0xPTw80bnJyMlNTUyPPsVAWBjOSIgoAAJy+Op1O7v7C3dl03qaBxq05vCZJcujeQ0sec+/BexfJ8aVsOv/ywXIc6dagQ199ZKBx9x74+kD7099IiujExEQO5+y89lnPbTLfjZ/8SFZPbOyb5VCO5rU/8Mzh5/jEtqyZmOib45t1f151zSVDz5Ekb/rUg3lCnywAALCYTedtyqs3v3ro87xx+xsXznH+5XnNM39+6DmS5A3bfrfJPGcC14gCAADQlCIKAABAU4ooAAAATSmiAAAANKWIAgAA0JQiCgAAQFOKKAAAAE0pogAAADS1atQBGC+zs7M5uLfmrZ881Ga+vTXnZbZvlgP7at718SNDz/HA3ppHa/8c+/clWz92dOg5kmT33uRQnywAAIynmZmZdDqdgcYc3396enrg+SYnJzM1NbUsOU4lS78ci1FEAQAATlGn08ndX/hiNl1w6ZLHrDnSPUH10Nf2DTTXvfvvXyTHzmxau2Ggx1zTO+Zy6Ot7lp5j366B5jiRIsp/Z2JiIgezJz/7rDVN5nvrJw/lvImJvln2lj35yR8c/rfpuz5+JBde1j/HmrI71//Q2UPPkXSPvK7vkwUAgPG16YJL85++/2VDn+f1f33zwjnWbshrr3nx0HPc+Kn3nvRY14gCAADQlCIKAABAU4ooAAAATSmiAAAANKWIAgAA0JQiCgAAQFOLFtFSyjtKKQ+WUu464b5fK6V8vZRyR+/f84cbEwAAgJViKUdE35nkunnuf0ut9erevw8tbywAAABWqlWL7VBr/UQp5crhRwEWMjs7m337kr+4rc18Dz+cHDs22zfL/n3Jp/+8Dj3H/oeS2aP9czy6P9n5wWNDz5Ekj+5JZg/PzTI7O5vsT469/0iTHNmdzB6a/zkBADgdnMo1oq8spdzZO3X3on47lVJuKKVsL6Vs37Vr1ylMBwAAwEqw6BHRPmaSvC5J7X18c5KfmW/HWutNSW5Kks2bNw//8AmsUBMTEznrrN354R9pM99f3JZceulE3yzHzt6d73teGXqOT/95zcTG/jm+uXp3nv6jbV53becHj2Viw9wsExMT2b1md8768ZNdUgdz7P1HMrF+/ucEAOB0cFK/vdVaH6i1Hq21Hkvy+0mesbyxAAAAWKlOqoiWUi474eaPJbmr374AAABwokXPIyul/HGSZydZX0r5WpL/kuTZpZSr0z019ytJ/t0QMwIAALCCLOVVc186z91vH0IWAAAAzgBtXuEDAAAAehRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYWfdVcAE4fs7Ozyf5Hc+QDO4c/2e5HM3totm+Ouv9AjmzdNvwcSeqeA5k9PH8WAGD8OCIKAABAU46IAqwgExMT2b3mG1n1oqcPfa4jH9iZifUTfXPsWX00q65/5tBzJMmRrdsysWH+LADA+HFEFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgqVWjDgAAwzI7O5u6f1+O3HJ7k/nqnr2ZPVwXzHL41g83yPFQZg8fXSDH/hy6devQc3Sz7M7s4UNN5gLg9OGIKAAAAE05IgrAijUxMZE9q0tWvfDZTeY7csvtmdhw2QJZzs7qF1w39ByHb/1wJjZsXCDHmqx5wfVDz5Ekh27dmokN65vMBcDpwxFRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaWjXqAADAmWd2djbH9h/IN259d5P5ju15MLOHH+ub5ej+A3n0lrcNPcfRPfdl9vDBoc8DMO4cEQUAAKApR0QBgOYmJiby0OqH88QXvKTJfN+49d2Z2HBR3ywPr96fc174iqHnePSWt2ViwwVDnwdg3DkiCgAAQFOKKAAAAE0pogAAADSliAIAANCUIgoAAEBTiigAAABNKaIAAAA0pYgCAADQ1KpRBwAAIJmdnc2R/Y9k/9Y3N5nvyJ6vZvbwufPmOLD/kWzb+oYmOQ7sviezh+bmAFY2R0QBAABoyhFRAIAxMDExkb2rv5ELrv/lJvPt3/rmTGx44rw5jq45lGde/5omObZtfUMm1q9pMhcwPhwRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFi2ipZR3lFIeLKXcdcJ9F5dSbiulfLn38aLhxgQAAGClWMoR0Xcmue5x9/1qko/VWp+a5GO92wAAALCoRYtorfUTSR563N3XJ7m59/nNSV60zLkAAABYoVad5LiNtdb7kqTWel8p5ZJlzAQAAJmdnc3+/Y/mllteN/S59uy5J4cPn9M3x4F9j+YPP3Lj0HMkyQMP3ZNHj87NMjs7m4N7H8nv3d4mx9f33pPzcu6822ZnZ/PIgUfyxu1vHHqOew7ck3NnF8pxMG/Y9rtDz9HN8vWcO3tek7lWuqG/WFEp5YZSyvZSyvZdu3YNezoAAADG3MkeEX2glHJZ72joZUke7LdjrfWmJDclyebNm+tJzgcAwBlmYmIiq1cfzgtf+J+HPtctt7wuGzas7ptj79mH81PPfe3QcyTJH37kxly4cW6WiYmJPJJD+blnt8nxe7ffmHMn1sy7bWJiIoeOHMqrN7966DneuP2NWbNQjqOP5DXP/Pmh50iSN2z73ayZmP/oLIM52SOityR5We/zlyXZujxxAAAAWOmW8vYtf5zk00m+vZTytVLKy5P8ZpIfKaV8OcmP9G4DAADAohY9NbfW+tI+m35ombMAAABwBhj6ixUBAADAiRRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYWfdXcYbln30O58ZMfGWjM/QcPJEkuPe/8gee6amLjQGMAAACWanZ2No/sP5DX//XNQ5/rnv3359zZR/rn2Lc/N37qvcPPsW9Xzi3fPKmxIymik5OTJzXuUOdgkmT1gKXyqomNJz0nAAAAy2skRXRqauqkxk1PTydJtmzZspxxAAAATsnExEQOHduX//T9Lxv6XK//65uzZmJt/xz1CXntNS8eeo4bP/XerJlYd1JjXSMKAABAU4ooAAAATSmiAAAANKWIAgAA0JQiCgAAQFOKKAAAAE0pogAAADSliAIAANDUqlEHGAf37juQGz+xbcn7P3Dw0STJxvPOGXieqyYW2n44b/rUgwM95oOPHEmSXHLuYF/Ke/cdzlMvn3/b7L6at37y0ECPt/tgTZKsP68MNG52X83TFnhOHthb866PH1ny4z3cy3HRgDke2Ftz4WX9t+/em2z92NGBHnPfwe7HtecNNCy79ybrF8gCAACnuzO+iE5OTg485lCnkyRZM/HkgcZdNdF/vpPJcWKWJ1w+2PinXj7/nCeb44FejvMmBhv/tGV+Tvb0clx42WBjL7xs+b82+3pZ1g+YZf0CWQAAYCU444vo1NTUwGOmp6eTJFu2bBlpjmFkGZccJ5tlXHIMKwsAAKwErhEFAACgKUUUAACAphRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKCpVaMOACzdww8nf3HbYGMOHOh+PP/8wee69NL+2/c/lHz6z+uSH++RXo5zB8yx/6EkG/tvf3RPsvODxwZ6zG/s63584trBsjy6J8mGPht3J8fef2SwB+zlyIA5sjvJ+oW2P5YjH9g5YJZv9rI8YYAcjy2Yo+45kCNbtw0Uo+57NElS1p4z2Lg9B/p/bQCAsaOIwmlicnLypMY98kgnSXLppYONv/TS/nOeTJbOwW6OKzYOOHbj8uZIks7+bpZv2zDg+A3zz3nSOfZ1c0yuH3D8+iE8J9/K8j+GX0xcAAASLklEQVSMNkfvazO54cmDDezztQEAxpMiCqeJqampkxo3PT2dJNmyZctIs4xLjmFkGZcc45RlXHIAAOPJNaIAAAA0pYgCAADQlCIKAABAU4ooAAAATSmiAAAANKWIAgAA0JQiCgAAQFOKKAAAAE2tGnUAABimumdvjtxy+2Bj9h1MkpS15w08VzZctsD2h3L41g8PmOVAL8v5A+R4KNmwcYHtu3Po1q0D5tjXy7F2sHF7dicb1s+77dieB/ONW9890OMd2/dwkuSstRcNNm7Pg8mGwcYAMDyKKAAr1uTk5EmN6+zvdMcvUCrnteGyvnOefJaDvSz9i+XcHBuHkGNfL8f8pbJ/lvXzznnyOR7q5RiwVG646KTnBGD5KaIArFhTU1MnNW56ejpJsmXLlhWXRQ4AxoFrRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgqVN6+5ZSyleSHEhyNMmRWuvm5QgFAADAyrUc7yP6g7XW3cvwOAAAAJwBnJoLAABAU6d6RLQm+WgppSb5P2utNy1DJgCApo7uuT+P3vK2gcYc27cnSXLW2nUDzZMNFyyw/avZv/XNA+U4uu/BJMnZay8ZbNyeryYbnjrvtgO77822rW8Y6PEe3fdAkuSctRsHGndg973J+qv6bt+z557ccsvrBnrMffvuT5KsXXvpksfs2XNPNmzon+OBh+7JH37kxoFyPHygm+Oi85ee4/hcF26cP8vX992b37t9sBy7D3a/NuvPG+xr8/V99+ZpE/2fk3sP3ps3bn/jQI/5wKPdLBvPWXqWew/em6vSPwenp1MtotfUWmdLKZckua2UsrPW+okTdyil3JDkhiTZtGnTKU4HALC8JicnT2pcZ/+u7vgFiuUcGy7oO9/J5zjcy/HEwQZueOq8c550jn2HkiRPXr9msIHrr1r252T//m6WDRtWL3nMhg3Ln2PPwW6OCzcuPUd3//mznGyO+zvdHOdODPa1edrE8j8nh3pZ1mxaepar0j8Hp69TKqK11tnexwdLKX+a5BlJPvG4fW5KclOSbN68uZ7KfAAAy21qauqkxk1PTydJtmzZIscQcoxTFjnGOwunp5O+RrSUcm4p5fzjnye5NsldyxUMAACAlelUjohuTPKnpZTjj/NHtdYPL0sqAAAAVqyTLqK11n9I8t3LmAUAAIAzgLdvAQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKCpU3n7FgAAgJG698DX84ZtvzvQmAce3Z0k2XjO+oHnuipP6799//15/V/fvPQcjzzUzXHuxYPl2H9/rsra/tv37cqNn3rvQI/5wCN7e1kuXHqOfbty1eXrBprnOEUUAAA4LU1OTp7UuEOd+5Mka5587kDjrsrT+s55MlkOdbqFeM0V/Uvl/DnWLmuObpZuEV0zQLG86vJ1Jz2fIgoAAJyWpqamTmrc9PR0kmTLli0jzTIuOYaVZSGuEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgKUUUAACAphRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgKUUUAACAphRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABoShEFAACgKUUUAACAphRRAAAAmlJEAQAAaEoRBQAAoClFFAAAgKYUUQAAAJpSRAEAAGhKEQUAAKApRRQAAICmFFEAAACaUkQBAABo6pSKaCnlulLKF0spd5dSfnW5QgEAALBynXQRLaWcneR3kzwvyXcmeWkp5TuXKxgAAAAr06kcEX1Gkrtrrf9Qaz2U5N1Jrl+eWAAAAKxUpdZ6cgNLeXGS62qtr+jd/ukkz6y1vrLfmM2bN9ft27cv+LgzMzPpdDrzbjt+/+Tk5LzbJycnMzU1tZT4S9Ivy7jkaJ1lXHIslGVccrTOMi45FsoyLjlaZxmXHKeSZVxyLHeWcclxKlnGJcdyZxmXHKeSZVxyLHeWcclxKlnGJcdyZxmXHKeSZVxyLHeWcclxKlkGyVFK+WytdfNi+61a0qP1mWOe++a02lLKDUluSJJNmzadwnTJk570pFMav1zGJUcyPlnkmGtcssgx17hkGZccyfhkkWOucckix1zjkkWOucYlixxzjUsWOeZqneVUjoh+X5Jfq7U+t3f71UlSa31jvzFLOSIKAADA6WmpR0RP5RrRv0ny1FLKt5VS1iR5SZJbTuHxAAAAOAOc9Km5tdYjpZRXJvlIkrOTvKPW+vllSwYAAMCKdCrXiKbW+qEkH1qmLAAAAJwBTuXUXAAAABiYIgoAAEBTiigAAABNKaIAAAA0pYgCAADQlCIKAABAU4ooAAAATSmiAAAANKWIAgAA0JQiCgAAQFOKKAAAAE0pogAAADSliAIAANCUIgoAAEBTiigAAABNKaIAAAA0pYgCAADQVKm1tpuslF1J7jnFh1mfZPcyxDlV45IjGZ8scsw1LlnkmGtcsoxLjmR8ssgx17hkkWOucckix1zjkkWOucYlixxzLUeWp9RaNyy2U9MiuhxKKdtrrZvl+CfjkkWOucYlixxzjUuWccmRjE8WOeYalyxyzDUuWeSYa1yyyDHXuGSRY66WWZyaCwAAQFOKKAAAAE2djkX0plEH6BmXHMn4ZJFjrnHJIsdc45JlXHIk45NFjrnGJYscc41LFjnmGpcscsw1LlnkmKtZltPuGlEAAABOb6fjEVEAAABOY2NbREspTy6lfLyUsqOU8vlSyi+csO1/K6V8sXf/fx1VllLKn5RS7uj9+0op5Y4R5bi6lPL/9XJsL6U8Y5g5Fsny3aWUT5dSPldKubWUcsGQczyxlPKZUsrf93L8eu/+byulbCulfLn3dVozohyvLKXcXUqppZT1w8ywhCzv6v3c3FVKeUcpZfWIcry9d9+dpZT3llLOG0WOE7b/t1LKwWFmWCxLKeWdpZR/PGE9uXpEOUop5fWllC/1frb//YhyfPKE52K2lPKBYeZYJMsPlVL+tpflr0opV40ox3N6Oe4qpdxcSlk1zBwn5Dm7lPJ3pZQP9m43XVsXyNF8bV0gS9O1dYEcTdfWfjlOuL/Z2tovS+u1dYEcTdfWBXI0X1sXyNJ0bV0gx6jW1q+U7u/Kd5RStvfuu7iUcltvfb2tlHLRiHL8RO//n2OllOG+em6tdSz/JbksyT/rfX5+ki8l+c4kP5jkL5I8obftklFledw+b07yv4/oOflokuf17n9+kttH+PX5myT/onf/zyR53ZBzlCTn9T5fnWRbku9N8p4kL+nd/9YkUyPK8T1JrkzylSTrh/11WSTL83vbSpI/HuFzcsEJ+/x2kl8dRY7e7c1J/u8kB0f8tXlnkhe3yLBIjn+T5P9KclZv21DX14W+Nifs874k/3qEz8mXknxH7/6fS/LOEeT4/iRfTfK03v2/keTljb5X/kOSP0rywd7tpmvrAjmar60LZGm6ti6Qo+na2i9H776ma+sCz0nTtXWBHE3X1oW+Nidsa7K2LvCcNF1b58uR7gG5Ua2tc9auJP/1+M9tkl9N8qYR5fiOJN+e5PYkm4c5/9geEa213ldr/dve5weS7EhyeZKpJL9Za/1mb9uDI8ySpPuXriT/a7r/AY0iR01y/Mjj2iSzw8yxSJZvT/KJ3m63Jflfhpyj1lqP/8V1de9fTfKcJO/t3X9zkheNIket9e9qrV8Z5twDZPlQb1tN8pkkV4wox/7kWz83T0r369U8Rynl7CS/leQ/DnP+pWRpNf8Sckwl+Y1a67HefkNdXxd7Pkop56f7szz0v9ovkKXp+tonx9Ek36y1fql3/9DX1iQppVyR5F8meVvvdknjtXW+HEkyirV1gSxN19YFcjRdW/vlGMXa2i/LKPTJ0XRtXSDH8W3N1tYFsjT/3XWeHOsygrV1Adenu64mjdbX+dRad9Rav9hirrEtoicqpVyZ7l8/tyV5WpJn9U4N+n9LKf98hFmOe1aSB2qtXx5Rjl9M8lullK8m2ZLk1a1yzJPlriQv7G36iSRPbjD/2aV7WvSD6S4inSR7a61Hert8LSf84aBVjlrrtsXGjCJL77Sxn07y4VHlKKX8QZL7kzw9yX8bUY5XJrml1nrfsOdfQpYkeX3vlLq3lFKeMKIck0n+Veme4v/npZSnjijHcT+W5GPHf8EeUZZXJPlQKeVr6f7c/GbrHOmWm9UnnCL14jRYW5P8Trpl4ljv9rqMYG2dJ8co9c3Scm3tl6P12tonx0jW1j5ZksZra58czdfWPjmOa7q29snSfG2dJ8fujGZtTbpF/KOllM+WUm7o3bfx+M9N7+MlI8rRzNgX0dK9xuF9SX6x9wOzKslF6Z4y9StJ3tP7698oshz30gz5aOgiOaaS/FKt9clJfinJ20eY5WeS/Hwp5bPpnrJ7aNgZaq1Ha61Xp/tX6Geke0rBnN1a5yilfNew5zzJLL+X5BO11k+OKket9d8kmUj3SPq/GkGOH0j3DyUtflFbLMt3pfvHo6cn+edJLk7yqhHleEKSb9RaNyf5/STvGFGO45qurX2y/FKS59dar0jyB+me8tg0R5L/MclLkryllPKZJAeSHFngIU5ZKeVHkzxYa/3siXfPF3cEOUZiCVmarK0L5Wi5ts6Xo5QykRGsrQs8J03X1gVyNF1bl/C92mxtXSBL07V1vhy9sxiarq0nuKbW+s+SPC/d35t/oNG8Y5VjrIto76+L70vyrlrr+3t3fy3J+3tnwXwm3b9qDP2FCvpkSe+i5h9P8ifDzrBAjpclOf75/5PuLy4jyVJr3VlrvbbW+j+nu8h1WmTpzb033fPZvzfJhSdccH5FGpzyMU+O61rN2c/js5RS/kuSDeleIzGyHL37jqb7c9PsNJgTcvxgkquS3F1K+UqSc0opd7fK8bgs19Xuqe61di85+IM0+hl+fI5019f39Tb9aZL/aUQ5UkpZl+7z8GetMsyT5XlJvvuEo7R/ku71mq1zXFdr/XSt9Vm11meke/nDsM/AuSbJC3s/H+9O9zS+30n7tXVOjlLKHw55zoGzNF5bF3xOGq6t832PfD6jWVvnfU5GsLb2+9q0XlsX+l5tvbbOl+XP0n5t7fc90nptTZLUWmd7Hx9M93viGUkeKKVcliS9jy0uP5wvRzNjW0R7RznfnmRHrfXEv5J8IN3FLqWUpyVZk+6h9VFkSZIfTrKz1vq1YWZYJMdskn/R+/w5afBD1C9LKeWS3sezkrw23RezGGaODaWUC3ufPyndr8eOJB9P9xSLpFvUt44gx85hzjlollLKK5I8N8lLj1+nMoIcXyy9V8brfQ+9IEN+nvrk+Gyt9dJa65W11iuTPFprHfor9i3wtTn+H09J95qQu0aRIyesr+muKV+a/xGGniPpHlX5YK31G8PMsEiWHUnW9v6vSZIf6d3XOsfOE9bWJ6R7VGeoa2ut9dW11it6Px8vSfKXtdafTOO1tU+OnxrmnINmab22zpcjyU+3Xlv7PB8XjWJtXeBr03RtXeD7tenausjPTdO1tc/36/VpvLYu8D3SdG3tzXVu6V6nm1LKuUmuTfd785Z019Wkze+u/XI00+Qlik/SNemeM/658k9vi/KadE9neEcp5a50T/t8We/QevMstdYPpfvN3OrUsX7Pyb9N8n/0/kr9jSQtzvHul+WppZSf791+f7p/gRymy5LcXLovjnBWkvfUWj9YSvlCkneXUm5M8ncZ/unK/XL8+3SvR7g0yZ2llA/VWl8xoixHktyT5NPd/5Pz/lrrb7TMke5fYD9Zum/rU5L8fbqnlg/TvM/HkOccKEsp5S9LKRvSfU7uSPKzI8rxV0neVUr5pSQH072Gp3mO3raXpM01QwtmKaX82yTvK6UcS/JwupcfjCLHb5XuqWVnJZmptf7lkHP086q0XVvnNaK1tZ+3pu3aOp+S7vdNy7X1dPCuxmtrP7+ZtmvrQlqvrXPUWo+MYG3t51dGsLZuTPKnvfViVZI/qrV+uJTyN+lecvjyJPem+0eDUeT4sXRPr9+Q5M9KKXfUWp87jABl+B0OAAAA/snYnpoLAADAyqSIAgAA0JQiCgAAQFOKKAAAAE0pogAAADSliAIAANCUIgoAAEBTiigAAABN/f/n3k98dz2S3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x648 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')\n",
    "qual_pos = defaultdict(list)\n",
    "for rec in recs:\n",
    "    for i, qual in enumerate(rec.letter_annotations['phred_quality']):\n",
    "        if i < 25 or qual == 40:\n",
    "            continue\n",
    "        pos = i + 1\n",
    "        qual_pos[pos].append(qual)\n",
    "vps = []\n",
    "poses = list(qual_pos.keys())\n",
    "poses.sort()\n",
    "for pos in poses:\n",
    "    vps.append(qual_pos[pos])\n",
    "fig, ax = plt.subplots(figsize=(16,9))\n",
    "sns.boxplot(data=vps, ax=ax)\n",
    "ax.set_xticklabels([str(x) for x in range(26, max(qual_pos.keys()) + 1)])\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# There is more..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do this to download the paired end data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Be careful as this will be 1GB of data (and fully optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-08-26 16:51:39--  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265_1.filt.fastq.gz\n",
      "           => ‘SRR003265_1.filt.fastq.gz’\n",
      "Resolving ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)... 193.62.192.8\n",
      "Connecting to ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)|193.62.192.8|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /vol1/ftp/phase3/data/NA18489/sequence_read ... done.\n",
      "==> SIZE SRR003265_1.filt.fastq.gz ... 502660639\n",
      "==> PASV ... done.    ==> RETR SRR003265_1.filt.fastq.gz ... done.\n",
      "Length: 502660639 (479M) (unauthoritative)\n",
      "\n",
      "SRR003265_1.filt.fa 100%[===================>] 479.37M  18.0MB/s    in 46s     \n",
      "\n",
      "2018-08-26 16:52:27 (10.4 MB/s) - ‘SRR003265_1.filt.fastq.gz’ saved [502660639]\n",
      "\n",
      "--2018-08-26 16:52:27--  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265_2.filt.fastq.gz\n",
      "           => ‘SRR003265_2.filt.fastq.gz’\n",
      "Resolving ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)... 193.62.192.8\n",
      "Connecting to ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)|193.62.192.8|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /vol1/ftp/phase3/data/NA18489/sequence_read ... done.\n",
      "==> SIZE SRR003265_2.filt.fastq.gz ... 484084218\n",
      "==> PASV ... done.    ==> RETR SRR003265_2.filt.fastq.gz ... done.\n",
      "Length: 484084218 (462M) (unauthoritative)\n",
      "\n",
      "SRR003265_2.filt.fa 100%[===================>] 461.66M  13.8MB/s    in 43s     \n",
      "\n",
      "2018-08-26 16:53:12 (10.7 MB/s) - ‘SRR003265_2.filt.fastq.gz’ saved [484084218]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#XXX change\n",
    "!rm -f SRR003265_1.filt.fastq.gz 2>/dev/null\n",
    "!rm -f SRR003265_2.filt.fastq.gz 2>/dev/null\n",
    "!wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265_1.filt.fastq.gz\n",
    "!wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265_2.filt.fastq.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of pairs: 9170808\n"
     ]
    }
   ],
   "source": [
    "f1 = gzip.open('SRR003265_1.filt.fastq.gz', 'rt', encoding='utf-8')\n",
    "f2 = gzip.open('SRR003265_2.filt.fastq.gz', 'rt', encoding='utf-8')\n",
    "recs1 = SeqIO.parse(f1, 'fastq')\n",
    "recs2 = SeqIO.parse(f2, 'fastq')\n",
    "cnt = 0\n",
    "for rec1 in recs1:\n",
    "    next(recs2)\n",
    "    cnt +=1\n",
    "print('Number of pairs: %d' % cnt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Only do the next cell on Python 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of pairs: 9170808\n"
     ]
    }
   ],
   "source": [
    "f1 = gzip.open('SRR003265_1.filt.fastq.gz', 'rt', encoding='utf8')\n",
    "f2 = gzip.open('SRR003265_2.filt.fastq.gz', 'rt', encoding='utf8')\n",
    "recs1 = SeqIO.parse(f1, 'fastq')\n",
    "recs2 = SeqIO.parse(f2, 'fastq')\n",
    "cnt = 0\n",
    "for rec1, rec2 in zip(recs1, recs2):\n",
    "    cnt +=1\n",
    "\n",
    "print('Number of pairs: %d' % cnt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
